1
$\begingroup$

I'm modeling stochastic chemical kinetics and the ItoProcess[] function has served me well. I am trying to write an efficient code to analyze many (Paths) trajectories of several different reagents (vector output)

If I have two chemical reactions:

A -> XL

A -> XD

I can write the Langevin equation for each species to simulate the reaction using a Weiner Process.

enter image description here

I am interested in the time course trajectories of A[t], XL[t], XD[t] as well as the terms XL[t]-XD[t] and (XL[t]-XD[t])/(XL[t]+XD[t]).

I can create an ItoProcess with two WeinerProcesses, and the variables above as the output (ai = initial particle numbers of a):

ai = 1000;    
proc = ItoProcess[
       {\[DifferentialD]A[t] == -2 A[t] \[DifferentialD]t - 
          A[t]^0.5` \[DifferentialD]w1[t] - 
          A[t]^0.5` \[DifferentialD]w2[t],
        \[DifferentialD]XL[t] == 
         A[t] \[DifferentialD]t + A[t]^0.5` \[DifferentialD]w1[t],
        \[DifferentialD]XD[t] == 
         A[t] \[DifferentialD]t + A[t]^0.5` \[DifferentialD]w2[t]},
       {A[t], XL[t], XD[t], XL[t] - XD[t],
        (XL[t] - XD[t])/(XL[t] + XD[t])},
       {{A, XL, XD}, {ai, 0, 0}},
       t,
       {w1 \[Distributed] WienerProcess[], 
        w2 \[Distributed] WienerProcess[]}];

I can see the output of one individual trajectory as a temporal function using:

tend = 200; (*end time *)
tint= 200/10^4; (*number of intervals*)
npaths = 1 (* number of trajectories *)   

td = RandomFunction[proc, {0.0001, tend, tint},npaths];
ListLinePlot[td]

enter image description here

This is all well and good for a scenario with one path, but I'd like to perform this analysis with thousands of paths. I would like to output each of the vector values A[t], XL[t], XD[t], XL[t]-XD[t] as an individual temporal data file (with multiple paths). That way I can use properties such as "SliceData"," SliceDistribution", "Paths", ect to analyze the data easily with a large number of paths.

I can pull out individual trajectories using a TableFunction like so:

xltable = Partition[Riffle[Table[i*tint, {i, 0, tend/tint}], td["Values"][[All, 1]]],2];
ListLinePlot[xltable]

but this becomes slow when npath is high and the ItoProcess is more complex.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.