How to do a simple parameter scan using Tellurium

Originally Posted on August 27, 2018 by hsauro

A common task in modeling is to see how a parameter influences a model’s dynamics. For example, consider a simple two reaction pathway:

-> S1 ->

where the first reaction has a fixed input of vo and the second reaction a first-order rate laws k1*S1. The task is to investigate how the time course of S1 is influenced by vo.

The script below defines the model, then changes vo in increments and plots the effect on the pathway via a time course simulation. To do the parameter scan we exploit plotArray. This initially prevents the plots from being shown using show=False. We make sure that each plot gets a different color using resetColorCycle=False, finally, we show the plot using show(). To make things more interesting we also add a legend entry for each plot.

Note we call reset each time we run a simulation to ensure that S1 is reset back to its initial condition.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

import tellurium as te import numpy as np   r = te.loada("""         J1: -> S1; vo;         S1 -> ; k1*S1         k1 = 0.5;         vo = 4; """)   for vo in np.arange (0, 5, 0.5):     r.reset()     r.vo = vo;     m = r.simulate(0, 10, 100)     te.plotArray(m, labels=['v0 = '+str(vo)], show=False, resetColorCycle=False) te.show()

The following figure shows the resulting plot:

Thanks to Kiri Choi for pointing out how to use plotArray in this way.

How to plot a grid of phase plots using Tellurium

Originally Posted on August 22, 2018 by hsauro

Let’s say we have a chemical network model where the species oscillate and you’d like to plot every combination of these on a grid. If so, then this code might be of help.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

import tellurium as te import roadrunner import pylab as plt   import tellurium as te r = te.loada (''' model feedback()    // Reactions:    J0: $X0 -> S1; (VM1 * (X0 - S1/Keq1))/(1 + X0 + S1 +   S4^h);    J1: S1 -> S2; (10 * S1 - 2 * S2) / (1 + S1 + S2);    J2: S2 -> S3; (10 * S2 - 2 * S3) / (1 + S2 + S3);    J3: S3 -> S4; (10 * S3 - 2 * S4) / (1 + S3 + S4);    J4: S4 -> $X1; (V4 * S4) / (KS4 + S4);     // Species initializations:   S1 = 0; S2 = 0; S3 = 0;   S4 = 0; X0 = 10; X1 = 0;     // Variable initialization:   VM1 = 10; Keq1 = 10; h = 10; V4 = 2.5; KS4 = 0.5; end''')   slist = r.getFloatingSpeciesIds() fig = plt.figure(figsize=(9,6)) fig.tight_layout() r.reset() m = r.simulate (0, 80, 2000, slist) count = 0 for i in range(4):     for j in range(4):         count += 1         if i >= j:             ax = fig.add_subplot(4, 4, count)             ax.set_xlabel (slist[i])             ax.set_ylabel (slist[j])             ax.plot (m[:,i], m[:,j])

This code will generate:

Here is an alternative that has removed the internal ticks and labels and arranges the column and row labels as contiguous. This is probably more than one might expect. This version plots all combinations including the transpose combinations.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52

import tellurium as te import roadrunner import pylab as plt   import tellurium as te r = te.loada (''' model feedback()    // Reactions:    J0: $X0 -> S1; (VM1 * (X0 - S1/Keq1))/(1 + X0 + S1 +   S4^h);    J1: S1 -> S2; (10 * S1 - 2 * S2) / (1 + S1 + S2);    J2: S2 -> S3; (10 * S2 - 2 * S3) / (1 + S2 + S3);    J3: S3 -> S4; (10 * S3 - 2 * S4) / (1 + S3 + S4);    J4: S4 -> $X1; (V4 * S4) / (KS4 + S4);     // Species initializations:   S1 = 0; S2 = 0; S3 = 0;   S4 = 0; X0 = 10; X1 = 0;     // Variable initialization:   VM1 = 10; Keq1 = 10; h = 10; V4 = 2.5; KS4 = 0.5; end''')        count = 0 slist = sorted (r.getFloatingSpeciesIds()) m = r.simulate (0, 80, 2000, slist) fig, axarr = plt.subplots(4,4,figsize=(9,6)) fig.subplots_adjust (wspace=0.0, hspace=0.0)   for i in range(4):     for j in range(4):         count += 1         ax = plt.subplot2grid ((4,4),(i,j))         if i==3:            ax.set_xlabel (slist[j])         else:            ax.set_xticklabels([])            ax.set_xticks([])         if j == 0:            ax.set_ylabel (slist[i])                       else:            ax.set_yticklabels([])            ax.set_yticks([])           ax.plot (m[:,i], m[:,j])                  # Comment the above line and uncomment these lines         # to only plot the lower half of the grid         #if count in [2,3,4,7,8,12]:         #   ax.set_visible(False)         #else:         #   ax.plot (m[:,i], m[:,j])

The following shows the steady-state oscillations. This was done by simulating the model twice and using the plotting the second set of simulation results.