Sunday, December 17, 2023

Animating pathway

Here is an example of runnign a simulation and plotting the results as an animation. The model is a 20 step linear pathway where at time zero, all concentrations are zero. The boundary species Xo is set to 10 and as the simualtion runs the species slowly fill up to their steady-state levels. We can capture the data and turn it intot an animation using the matplotlib.animation from matlpotlib. At the end we convert the frames into a mp4 file. I generated the model using hte buildnetworks class in teUtils

import tellurium as te
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, FFMpegFileWriter

r = te.loada('''
J1: $Xo -> S1; e1*(k10*Xo - k11*S1); 
J2: S1 -> S2; e2*(k20*S1 - k21*S2); 
J3: S2 -> S3; e3*(k30*S2 - k31*S3); 
J4: S3 -> S4; e4*(k40*S3 - k41*S4); 
J5: S4 -> S5; e5*(k50*S4 - k51*S5); 
J6: S5 -> S6; e6*(k60*S5 - k61*S6); 
J7: S6 -> S7; e7*(k70*S6 - k71*S7); 
J8: S7 -> S8; e8*(k80*S7 - k81*S8); 
J9: S8 -> S9; e9*(k90*S8 - k91*S9); 
J10: S9 -> S10; e10*(k100*S9 - k101*S10); 
J11: S10 -> S11; e11*(k110*S10 - k111*S11); 
J12: S11 -> S12; e12*(k120*S11 - k121*S12); 
J13: S12 -> S13; e13*(k130*S12 - k131*S13); 
J14: S13 -> S14; e14*(k140*S13 - k141*S14); 
J15: S14 -> S15; e15*(k150*S14 - k151*S15); 
J16: S15 -> S16; e16*(k160*S15 - k161*S16); 
J17: S16 -> S17; e17*(k170*S16 - k171*S17); 
J18: S17 -> S18; e18*(k180*S17 - k181*S18); 
J19: S18 -> S19; e19*(k190*S18 - k191*S19); 
J20: S19 -> $X1; e20*(k200*S19 - k201*X1); 

k10 = 0.86;  k11 = 0.24
e1= 1; k20 = 2.24;  k21 = 0.51
e2= 1; k30 = 2.00;  k31 = 0.49
e3= 1; k40 = 2.54;  k41 = 0.83
e4= 1; k50 = 3.05;  k51 = 0.97
e5= 1; k60 = 1.29;  k61 = 0.10
e6= 1; k70 = 1.94;  k71 = 0.20
e7= 1; k80 = 3.01;  k81 = 0.17
e8= 1; k90 = 0.64;  k91 = 0.19
e9= 1; k100 = 1.96;  k101 = 0.42
e10= 1; k110 = 2.47;  k111 = 0.63
e11= 1; k120 = 4.95;  k121 = 0.51
e12= 1; k130 = 4.78;  k131 = 0.70
e13= 1; k140 = 3.77;  k141 = 0.78
e14= 1; k150 = 1.93;  k151 = 0.37
e15= 1; k160 = 3.87;  k161 = 0.96
e16= 1; k170 = 0.83;  k171 = 0.37
e17= 1; k180 = 1.82;  k181 = 0.69
e18= 1; k190 = 2.62;  k191 = 0.83
e19= 1; k200 = 4.85;  k201 = 0.67
e20= 1; Xo = 10.00
X1 = 0
S1 = 0; S2 = 0; S3 = 0; S4 = 0; 
S5 = 0; S6 = 0; S7 = 0; S8 = 0; 
S9 = 0; S10 = 0; S11 = 0; S12 = 0; 
S13 = 0; S14 = 0; S15 = 0; S16 = 0; 
S17 = 0; S18 = 0; S19 = 0; 
''')

# I'm running this on windows so I had to install the ffmpeg binary
# in order to save the video as an mp4 file
# for more info go tot his page:
# https://suryadayn.medium.com/error-requested-moviewriter-ffmpeg-not-available-easy-fix-9d1890a487d3
plt.rcParams['animation.ffmpeg_path'] ="C:\\ffmpeg\\bin\\ffmpeg.exe"

endTime = 18
fig, ax = plt.subplots(figsize=(8, 5))
ax.set(xlim=(-1, r.getNumIndFloatingSpecies()), ylim=(0, 18))

m = r.simulate(0, endTime, 100)

label = ax.text(15.6, 19, 'Time=', ha='center', va='center', fontsize=18, color="Black")
bars = ax.bar(r.getFloatingSpeciesIds(), m[0,1:], color='b', alpha = 0.5)
    
  def barAnimate(i):
    label.set_text('Time = ' + f'{m[i,0]:.2f}')

    for bar, h in zip(bars, m[i,1:]):
        bar.set_color ('r')
        bar.set_height(h)

anim = FuncAnimation(fig, barAnimate, interval=50, frames=100, repeat=False)

# You can save themp4 file anywhere you want
anim.save('c:\\tmp\\animation.mp4', fps=30)

plt.draw()
plt.show()

This is the resulting video:

Wednesday, December 13, 2023

Update to teUtils

 I just updated teUtils to version 2.9

UPDATE: The documentation was broken, now fixed. This is a set of utilities that can be used with our simulation environment Tellurium

The update includes some updates to the build synthetic networks functionality. 

The new version adds 'ei*(' terms to mass actions rate laws, eg
S1 -> S2; e1*(k1*S1 - k2*S2)

This makes it easier to compute control coefficients with repect to 'e' I also added a new mass-action rate of the form:

v = k1*A(1 - (B/A)/Keq1)

These can be generated using the call:

model = teUtils.buildNetworks.getLinearChain(10, rateLawType="ModifiedMassAction")

Useful if you want to more easily control the equilbrium constant for a reaction. Here is an example of a three step random linear chain:

model = teUtils.buildNetworks.getLinearChain(10, rateLawType="ModifiedMassAction")
print (model)

J1: $Xo -> S1; k1*Xo*(1 - (S1/Xo)/Keq1); 
J2: S1 -> S2; k2*S1*(1 - (S2/S1)/Keq2); 
J3: S2 -> $X1; k3*S2*(1 - (X1/S2)/Keq3); 

k1 = 1.42;  Keq1 = 3.61
e1 = 1; k2 = 1.42;  Keq2 = 7.91
e2 = 1; k3 = 3.50;  Keq3 = 5.64
e3 = 1; Xo = 5.00
X1 = 0
S1 = 1E-6; S2 = 1E-6; 

Monday, December 4, 2023

Executing tellurium models in a sphinx document

I needed to be able to generate simulation plots for a sphinx document but I didn't want to have to generate the plots separately and then include them manually.  I wanted the simulations done from within sphinx so that they would be automatically included when building the document. The key to this is to add the following to the list of extensions in the sphinx conf.py file:

extensions.append ('matplotlib.sphinxext.plot_directive')

To use it in a sphinx document use the plot directive, for example:

.. plot::
   :include-source:

   import tellurium as te
   r = te.loada ('''A -> B; k1*A; k1=0.1; A = 10''')
   m = r.simulate (0, 40, 100)
   r.plot()

It assumes you have matplotlib installed in your python setup (I am using Windows 10, Python 3.11). Further  information can be found here:

https://matplotlib.org/stable/api/sphinxext_plot_directive_api.html



How to execute python code in a sphinx document

This is a note for myself but others might find useful.  I needed to execute python in a sphinx document at build time. 

The package that comes up in a search doesn't work with Python 3.x. 

Instead I found 


which does work with Python 3.11. (I'm running Windows 10)

Install it by typing the following on the command line

pip install sphinx-execute-code-python3 

To use it don't forget to add this line to your sphinx conf.py file:

extensions.append('sphinx_execute_code')

Here is an example in a rst file:

.. execute_code::

    a = [1,2,3]
    b = [4,5,6]
    print ('Printing a list: ', a + b)

The docs have various options you can use.