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}},
       {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];

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];

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



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.