Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Has anyone figured out how to coerce ContinuousMarkovProcess to work with a state vector, transition matrix, and rate (propensity) vector, and RandomFunction to return pairs of state vector and time, so one can model chemical kinetics? The state vector contains the number of molecules of each species, and there are an infinite number of states in the traditional sense, so the state graph is infinite.

I have code for version based on version 8, but would like to make use of the new functions in version 9.

Here is a simple example of Lotka-Voltera system:

LT reactions

ODES:

{Derivative[1][Subscript[Y, 1]][t] == 
  Subscript[c, 1] Subscript[Y, 1][t] - 
   Subscript[c, 2] Subscript[Y, 1][t] Subscript[Y, 2][t], 
 Derivative[1][Subscript[Y, 2]][
   t] == -Subscript[c, 3] Subscript[Y, 2][t] + 
   Subscript[c, 2] Subscript[Y, 1][t] Subscript[Y, 2][t]}

Initial conditions:

{Subscript[Y, 1][0] == 4, Subscript[Y, 2][0] == 10}

Parameters:

{Subscript[c, 1] -> 1., Subscript[c, 2] -> 0.1, 
 Subscript[c, 3] -> 0.1}

Transition Matrix:

{{1, 0}, {-1, 1}, {0, -1}}

Propensities:

{Subscript[c, 1] Subscript[Y, 1][t], 
 Subscript[c, 2] Subscript[Y, 1][t] Subscript[Y, 2][t], 
 Subscript[c, 3] Subscript[Y, 2][t]}

Initial state vector:

{4, 10}

Thanks,

Bob

share|improve this question
    
Bob, a belated welcome to MMA.SE! –  Yves Klett Jan 30 at 15:46
    
Bob, meet Bob –  shrx Jan 30 at 15:57
    
Partial answer: If the propensities were constant, then one can used ContinuousMarkovProcess to find the series of reactions in the chain, then use the reactions and the stoich. matrix to construct the series of state vectors, then use Table[RandomReal[ExponentialDistribution[totalPropensity]], {Length[chain]}] to get the times. However, there are two issues. –  Bob Nachbar Feb 1 at 13:47
    
The first, and most important, is that the initial assumption above is false! The propensities depend on the state vector. Second, two different random sequences are used ({state1, state2, ...} and {time1, time2, ...} instead of a single interleaved sequence ({time1, state1, time2, state2, ...}). –  Bob Nachbar Feb 1 at 13:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.