Penultimate Version of Iridium Desktop Released. 

I believe I have finally reached the penultimate version of the Iridium Desktop simulator for Windows and Mac OS. The next version will be 1.0

The Iridium Desktop App is designed to be a productive and easy-to-use interactive SBML-compatible simulation environment for biochemical models. All code is open source  under the MIT license.

Mac OS and Windows binaries can be obtained from 

https://github.com/sys-bio/IridiumSimulator/releases

Click on assets to either download the Mac OS dmg file or a zip file for Windows. To use teh zip file, unzip the contents anywhere on your hard drive and run the Iridium executable. 

There is also a non-install Web Version at:

https://sys-bio.github.io/WebIridium/

Here are some screenshots of the Desktop version.

The main window is split into three sections, left to right: Control Panel, Input Panel, Output Panel. The control panel is switchable between different analysis types; these currently include Time Course simulation, Parameter Scanning, and Steady State. A variety of built-in models come with the distribution (Accessible in the bottom-left drop-down combo)

Time Course Simulation

The screenshot below shows the app operating in Time Course simulation mode. The central input panel has the antimony/SBML editor at the top and the optional slider panel at the bottom. Simulation results are shown on the plot panel to the right. 

Simulations are live in the sense that movement in any slider results in an instant change to the plotting output. All output is rendered in high definition. Output can also be captured as raw data by selecting the Text View tab on the output panel

Parameter Scan Panel

The following screenshot shows the parameter scan panel operation. Parameter scans can be with respect to any parameter of the model. The output measure can be any observable, but can also be limited to time course overlays, peak values, endpoints, and time to peak.  The slider panel is also active during a parameter scan.

Steady State Panel

The steady state panel is used to compute the steady state panel of the model. It can display the steady-state concentrations, Jacobian, and sensitivities of metabolic control analysis. The slider panel is also active in this mode. The steady state panel can display grid data or 3D plot.

Import/Export

SBML models (such as BioModels) can be imported or exported. When imported, they are automatically converted to Antimony form for easy editing.

Reproducibility

Both Time Course and Parameter Scan panels can export the last analysis as a Python script. This allows any parameter scan or time course simulation to be reproducible separately from the Iridium app.

Two sister applications are also under development:

1. Bifurcation Analysis Tool (Windows edition complete)

https://github.com/hsauro/BifurcationTool

2. A Model Fitting App for Windows and Mac OS

This is currently under development but should be released before October 2026

1. Fit SBML models to data using Antimony to edit models

2. Supports multiple data sets as well as time-course and dose-response data.

3. Implements 4 optimizers: Levenberg-Marquardt, Nelder-Mead, Swarm Particle, Differential Evolution.

3. Implements standard linear error uncertainty, profile likelihood, and bootstrap methods for uncertainty estimation.

4. Plots of the fitted model, residuals, profile likelihoods, and bootstrap estimations.

5. Differential Evolution and Particle Swarm can use multiple cores on a computer for higher performance.

6. Import/Export of PEtab.

All work on this page was supported by the NIH Center for Reproducible Biomedical Modeling

I created a number of cheat sheets for a class I did this quarter. The cheat sheets cover the topic of metabolic control. They include: 

1. General Metabolic Control 
2. Advanced Metabolic Control 
3. Linear Chains 
4. Branch Points 
5. Moiety Conserved Cycles

They are available a this GitHub repo:

https://github.com/hsauro/CheatSheets

Both the original LaTeX and pdfs files are provided.

The absolutely awful layout and render standard for SBML.

This summer I've spent 3 or 4 weeks trying to come to grips with the layout/render standard for SBML. This allows a modeler to specify a diagram, such as a metabolic pathway, alongside the model. It was proposed quite a few years ago but very few if any implementations have been made. I now understand why. The standard tries to do everything and as result it’s virtually impossible to fully implement, resulting in poor reproducibility (which was the whole point). There is one sublevel which is the layout layer but that gives very little in return. Access to the standard is via libsbml but the API is incredibly hard to use. It took me 250 lines of python to specify two reactions and three species. And even then I still hadn’t specified any colors, or thicknesses etc. To fully specify a two-reaction map might require up to 500 lines of python code. In C or C++ perhaps even more. One of my students has written a high-level API but even that is very hard to use and I gave up on it. Neither Claude nor ChatGPT can create working code using this ‘high-level’ API. Lucian Smith has written a layout/render extension to his Antimony language but even that has only partial support, but it is easier to use. There is some tooling out there but some of it is incomplete and the rest is hard to use. Here are some points worth making:

Poor API Design: The libSBML API requires creating multiple objects (BoundingBox, Point, Dimensions, Curve, LineSegment) just to draw a simple line. For example, to draw an arrow it should be sufficient to call draw_arrow(start, end).

Verbose XML Output: the standard is specified using XML which is slowly becoming an archaic format. While the model portion of an SBML model is readable (other than the MathML), the layout and render extension might as well be written in Sanskrit.

Tool Inconsistencies: Different viewers interpret the same SBML layout differently, making it unreliable for consistent visualization and reproducibility.

Missing Abstractions: There's no high-level concept of "draw pathway from A to B" - you have to manually construct every geometric primitive. The irony is that biochemical pathways are conceptually simple (nodes connected by directed edges), but the SBML standard makes them extraordinarily complex to represent.

The question is where next? SBML was developed almost 25 years ago. At the time we chose XML as the format carrier, and it was a good idea but today we have easier to manage formats, with better software support such as YAML, JSON etc. If we were to create SBML today it probably be something other than XML. The use of SBML to specify the visualization component was, however, a bad decision, or at least the specificaion is bad. Yes, SVG does it, but how long has it taken for SVG to become widely available? Even Google couldn't render SVG 1.1 until 2008. Even established vector drawing apps still don’t fully support it. Mobile support is also still spotty. With full industry backing it has taken a long time for SVG to become more mainstream.

In contrast to commercial settings, academic software development is heavily resource constrained, and the authors of the layout/render were perhaps a little optimistic that we'd be able to implement something as complex as SBML layout/render.

Layout and Render first came out in 2006 and the fact that there is hardly any support for it tells us the standard was too difficult to implement. The SBGN community shied away from it, probably because it was too complex. One of the best-practice rules we tried to develop during the development of SBML was that alongside a proposed standard there had to be at least one implementation that could exercise the standard to make sure it was a practical proposition. This happened with the model portion of SBML and showed us that software could be written without too much effort. The same applied I believe to the FBC extension. However, I don’t recall the same happening with the layout and render extension and this might explain the lack of implementations. Interestingly, the SBGN community didn’t even want to use it and instead developed their own ML.

So where do we go from here? Is it time to propose a successor to SBML that is easy to read and write and can incorporate extensions that can be implemented by the academic community?

For those interested, here is some example python code that tries unsuccessfully to create two reaction arcs (Yes each reaction, eve a ui-uni, has to have a minium of two curves). Ignore the silly if statements at the start, this was code under construction and no yet finalized but I gave up in the end. Note, this code just creates one reaction. Species and text creation is just as verbose.

def add_connections_to_reaction_glyph(reaction_glyph, reaction_id, species_positions, reaction_x, reaction_y):
    """
    Add reaction connections with proper Point syntax.
    """
    print(f"Adding connections for {reaction_id}...")
    
    # Define connections
    if reaction_id == 'J1':
        reactant_id = 'S1'
        product_id = 'S2'
    elif reaction_id == 'J2':
        reactant_id = 'S2'
        product_id = 'S3'
    else:
        return
    
    # Get species positions
    reactant_x, reactant_y, reactant_w, reactant_h = species_positions[reactant_id]
    product_x, product_y, product_w, product_h = species_positions[product_id]
    
    # Create reactant connection (species -> reaction)
    reactant_ref = reaction_glyph.createSpeciesReferenceGlyph()
    reactant_ref.setId(f"{reaction_id}_{reactant_id}_reactant")
    reactant_ref.setSpeciesGlyphId(f"{reactant_id}_glyph")
    reactant_ref.setRole(libsbml.SPECIES_ROLE_SUBSTRATE)
    
    # Create curve with proper Point syntax
    curve1 = libsbml.Curve()
    line_segment1 = libsbml.LineSegment()
    
    # Start point: right edge of reactant species
    start_point1 = libsbml.Point()
    start_point1.setX(reactant_x + reactant_w)  # Right edge
    start_point1.setY(reactant_y + reactant_h/2)  # Center height
    start_point1.setZ(0)
    
    # End point: reaction center
    end_point1 = libsbml.Point()
    end_point1.setX(reaction_x + 5)  # Center of 10x10 reaction
    end_point1.setY(reaction_y + 5)
    end_point1.setZ(0)
    
    line_segment1.setStart(start_point1)
    line_segment1.setEnd(end_point1)
    curve1.addCurveSegment(line_segment1)
    reactant_ref.setCurve(curve1)
    
    # Create product connection (reaction -> species)
    product_ref = reaction_glyph.createSpeciesReferenceGlyph()
    product_ref.setId(f"{reaction_id}_{product_id}_product")
    product_ref.setSpeciesGlyphId(f"{product_id}_glyph")
    product_ref.setRole(libsbml.SPECIES_ROLE_PRODUCT)
    
    # Create curve
    curve2 = libsbml.Curve()
    line_segment2 = libsbml.LineSegment()
    
    # Start point: reaction center
    start_point2 = libsbml.Point()
    start_point2.setX(reaction_x + 5)
    start_point2.setY(reaction_y + 5)
    start_point2.setZ(0)
    
    # End point: left edge of product species
    end_point2 = libsbml.Point()
    end_point2.setX(product_x)  # Left edge
    end_point2.setY(product_y + product_h/2)  # Center height
    end_point2.setZ(0)
    
    line_segment2.setStart(start_point2)
    line_segment2.setEnd(end_point2)
    curve2.addCurveSegment(line_segment2)
    product_ref.setCurve(curve2)

What were they thinking?

This paper popped into one of my feeds today:

Breakdown and repair of metabolism in the aging brain

I think the paper should have been titled "How not to publish a large model" The paper publishes a model is large but the way they deploy to the commuity is insane.

To save you hunting for the model, this link is to the GitHub repo

In summary the paper describes a kinetic model of brain metabolism with:

183 processes, which include:

95 enzymatic reactions

19 transport processes (across cell and mitochondrial membranes)

69 other processes (related to ionic currents, blood flow, and other non-enzymatic processes)

Additionally: The model uses 151 differential equations to simulate the dynamics of molecular concentrations.

So its large, but what's really a problem is the model is essentially inaccessible. The entire model is built using one huge Julia program. All the biology has been subsumed into a large set of difficult to read differential equations. There is no sharable SBML model so this won't go to Biomodels and reusing it will be very difficult. Why is this a problem? It means other researchers cannot build on what was undoubtedly, a huge amount of work. I took a screen shot of a small fragment of the Julia program so you can see what you're up against:

Another way to find unstable steady states

James Glazier recently told me of a trick he uses to find unstable steady states. Consider the following model which has two stable and one unstable steady state (it's a bistable system using postiive feedback). Note that $Xo means the species Xo is fixed.

import tellurium as te

r = te.loada('''
    $Xo -> S1; (0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')

If we run a simulation of this system it evolves to one of the stable steady states, in this case 0.144635. If we set the initial conditon to S1 = 10, we can also get the other stable state at S1 = 1.3095. There is the code to do that:

import tellurium as te

r = te.loada('''
    $Xo -> S1; ( 0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')

r.steadyState()
print (r.S1)

# Find theother steady state
r.S1 = 10
r.steadyState()
print (r.S1)

But how can we find the unstable one? There as an old trick where if one integrates backwards in time, stable states became repelers and unstable states attractors. However we don't allow someone to specfiy a start time that is bigger then the eed time. Instead James Glazier realized one chould just a put minus sign in front of every rate law to mimic the same effect. For example, a simulation of the following modiified model:

import tellurium as te

r = te.loada('''
    $Xo -> S1; -(0.1 + k1*S1^4/(k2+S1^4));
    S1 ->; -k3*S1;
    
    k1 = 0.9; k2 = 0.3; k3 = 0.7; S1 = 0.5;
''')
m = r.simulate (0, 50, 100)
print (r.S1)

will yield the unstable state at S1 = 0.68256.

Generating random networks

Here is a simple script that will generate a large number of random mass-action models and plot simulations of each in a grid.

Each grid has a number written on it so that you can grab the associated model if you find some interesting beahvior. In this case we see model 81 is oscillatory. You'll need the teUtils package for this in order to access the random network generator.

import tellurium as te
import roadrunner
import teUtils as tu
import matplotlib.pyplot as plt
import numpy as np

numRows = 15
numCols = 15
plt.subplots(numRows, numCols, figsize=(19,16)) 
count = 0
models = []
while count < numRows*numCols:
    model = tu.buildNetworks.getRandomNetwork(10,20)
    r = te.loada(model)
    try:
        m = r.simulate (0, 160, 200)
        try:            
            models.append (model)
            if count % 20 == 0:
               print (count)
            ax = plt.subplot (numRows, numCols, count+1)
            ax.tick_params(axis='both', which='major', labelsize=7)
            ax.tick_params(axis='both', which='minor', labelsize=7)
            ax.set_xlabel('Time')    
            te.plotArray(m, show=False)
            maxy = ax.get_ylim ()
            ax.text (50, maxy[1]/2, str(count), fontsize=14)            
            count = count + 1
        except:
            # failed to find a steady state so probably a bad model
            pass
    except:
        print ('Something very wrong with the model')
plt.show()

There is one run I did and you'll model 81 ihas some intersting dynamics:

This is model 81 pulled out and resimulated to show the dynamics more clearly:

import tellurium as te

# Get the model using:
#   print (models[81])

# Then copy and paste the model as below:

r = te.loada("""
var S0, S1, S2, S4, S5, S6, S7, S8, S9
ext S3;
J0: S9 + S8 -> S6; E0*(k0*S9*S8);
J1: S6 -> S7 + S9; E1*(k1*S6);
J2: S2 + S0 -> S8; E2*(k2*S2*S0);
J3: S4 + S4 -> S6; E3*(k3*S4*S4);
J4: S5 -> S2 + S1; E4*(k4*S5);
J5: S9 -> S4 + S0; E5*(k5*S9);
J6: S1 -> S0 + S5; E6*(k6*S1);
J7: S0 -> S4 + S4; E7*(k7*S0);
J8: S3 -> S9; E8*(k8*S3);
J9: S5 -> S1 + S9; E9*(k9*S5);
J10: S5 -> S0 + S4; E10*(k10*S5);
J11: S5 -> S4 + S9; E11*(k11*S5);
J12: S5 + S1 -> S0; E12*(k12*S5*S1);
J13: S2 + S9 -> S4; E13*(k13*S2*S9);
J14: S1 -> S8; E14*(k14*S1);
J15: S8 + S0 -> S1; E15*(k15*S8*S0);
J16: S7 -> S8; E16*(k16*S7);
J17: S8 + S6 -> S7; E17*(k17*S8*S6);
J18: S4 -> S6 + S6; E18*(k18*S4);
J19: S5 -> S1 + S4; E19*(k19*S5);

k0 = 0.6158; k1 = 0.0524;  k2 = 0.7206;  k3 = 0.0261
k4 = 0.4946; k5 = 0.2428;  k6 = 0.3249;  k7 = 0.4854
k8 = 0.6743; k9 = 0.6320;  k10 = 0.8954; k11 = 0.435
k12 = 0.580; k13 = 0.0298; k14 = 0.850;  k15 = 0.342
k16 = 0.556; k17 = 0.3221; k18 = 0.584;  k19 = 0.681

E0 = 1;  E1 = 1;  E2 = 1;  E3 = 1
E4 = 1;  E5 = 1;  E6 = 1;  E7 = 1
E8 = 1;  E9 = 1;  E10 = 1; E11 = 1
E12 = 1; E13 = 1; E14 = 1; E15 = 1
E16 = 1; E17 = 1; E18 = 1; E19 = 1

S3 = 1
S0 = 6; S1 = 3; S2 = 3; S4 = 2
S5 = 6; S6 = 2; S7 = 6; S8 = 2
S9 = 1
""")

m = r.simulate (0, 180, 100)
r.plot()

Here is the simulation from running the above code. Note that not every reaction is contributing to this behavior. You can remove J2, J3, J4, J6, J9, J10, J11, J12, J13, J14, J15, and J19 and the system will still oscillate.

What you can now is use another function in teUtils to look at a single model and plot a grid of simulations using random parameter values for that model. To do this just call plotRandSimGrid as follows. We apply it to model 81:

r = te.loada (models[81])
tu.plotting.plotRandSimGrid (r, endTime-500, ngrid=8, maxRange=1)

maxRange limits the range of random parmateer values. In this model, if the range is bigger, sometimes we get bad parameter sets which don't simulate. A maxRange of 1 seems to gives simulatable models. Most of the time the defaults are sufficient and you only have to pass in the roadrunner object, r. Here is a sample run:

Relationship of fluxes to enzyme levels in a metabolic pathway

To get hold of me it is best to contact me via my email address at: hsauro@uw.edu

Someone asked me the other day what the relationship was between the steady-state fluxe through a reaction and the coresponding level of enzyme. Someone else suggested that there would be a linear, or proportional relatinship between a flux and the enzyme level. However, this can’t be true, at least at steady. Considder a 10 step linear pathway. At steady-state each step in the pathway will, by defintion, carry the same flux. This is true even if each step has a different enzyme level. Hence the relationship is so simple. In fact the flux a given step carries is a systemic properties, dependent on all steps in the pathway. As an experiment I decided to do a simulation on some synthetic netowrks with random parameters and enzyme levels. For this exmaple I just used a simple rate law of the form: $$ v = e_i (k_1 A - k_2 B) $$ For a bibi reaction, A + B -> C + D, the coresponding rate law would be: $$ v = e_i (k_1 A B - k_2 C D) $$ a similar picture would be seen for the unibi and biuni reactions. Using our teUtils package I generated random networks with 60 species and 150 reactions. The reactions allowed are uiui-uinbi, biui or bibi. I then randomized the values for the enzymne levels $e_i$ and computed the steady-state flux. I used the following code to do the analysis. I have a small loop that generates 5 random models but obviously this number can be changed. I generate a random model, load the model into roadrunner, randomize the values for the $e_i$ parameters between 0 and 10, compute the steady-state (I do a presimulation to help things along) and collect the corresponding $e_i$ and flux values. Finally I plot each pair in a scatter plot.

import tellurium as te
import roadrunner
import teUtils as tu
import matplotlib.pyplot as plt
import random


for i in range (5):
    try:
      J = []; E = []
      antStr = tu.buildNetworks.getRandomNetwork(60, 150, isReversible=True)
      r = te.loada(antStr)

      n = r.getNumReactions()
      for i in range (n):      
          r.setValue ('E' + str (i), random.random()*10)    
      m = r.simulate(0, 200, 300)
      r.steadyState()

      for i in range (n):
          J.append (abs (r.getValue ('J' + str (i)))) 
          E.append (r.getValue ('E' + str (i)))        
      
      plt.figure(figsize=(12, 8))
      plt.plot(E, J, '.')
        
    except:
        print ('Error: bad model')

The results for five random networks is shown below. Note the x axis is the enzyme level and the y axis the corresponding steady-state flux through that enzyme. It's intersting to see that there is a rough correlation between enzyme amount and the corresponding flux, but its not very strong. Many of the points are just scattered randomly with some showing a definite correlation. The short answer is the realtinship is not so simple.

Experimenting with foreach loops in TikZ to draw biochemical pathways

Here is an example of a foreach loop in TikZ that can be used to draw an arbitrary long linear chain of reactions.

By setting the value of N in the following TikZ code, you can get any length linear chain. Obviously, you're limited by the width of the page. It uses the xifthen package to provide a conditional that is used to print out the last species, which is X1. There could be a better way of doing this, but this works. For example, use newcommand instead of def

 

\documentclass{article}

\usepackage{tikz}
\usetikzlibrary{arrows}
\usetikzlibrary{calc}
\usepackage{xifthen}

\begin{document}


\begin{tikzpicture}[>=latex', node distance=2cm]     
\node (X0) {$X_o$};
\pgfmathsetmacro{\N}{4}
\foreach \x in {0,...,\N}
{
  \pgfmathtruncatemacro{\nextval}{\x+1}
  \ifthenelse{\x = \N}
        {\def\speciesName{$X_1$}}
        {\def\speciesName{\large $x_\nextval$}}
  
  \node [right of = X\x] (X\nextval) {\speciesName};   
  \draw [->,ultra thick,blue] (X\x) -- node[above, black] {$v_{\nextval}$} (X\nextval);
}
\end{tikzpicture}

\end{document}

Here are some examples for N = 0, N = 2 and N = 4

MCA Rediscovered

It looks like someone has rediscovered metabolic control analysis (MCA).

A structural approach to understanding enzymatic regulation of chemical reaction networks Biochem J (2022) 479 (11): 1265–1283. Atsushi Mochizuki

The analysis is exactly the same as MCA but uses different symbols and they focus on the unscaled sensitivities instead. eg the r symbols are the unscaled elasticities. The core equation (4) can be found in equation (1) of the appendix of the following paper, and I am sure its been published elsewhere too:

Taking enzyme kinetics out of control; putting control into regulation. Jan-Hendrik S. HOFMEYR, Athel CORNISH-BOWDEN and Johann M. ROHWER, Eur. J. Biochem. 212, 833-837 (1993).

However, unlike the original MCA, the latest reincarnation doesn't include support for conserved moieties so as it stands it's somewhat limited. Note equation (1) in the above paper includes additional terms to take into account any conserved moieties. 

What is more concerning is that the reviewers of the paper never spotted this duplication of work. 

Some textbooks on MCA: 

Systems Biology: An Introduction to Metabolic Control Analysis (2018)

Understanding the Control of Metabolism (Frontiers in Metabolism Ser. ; No. 2)) 1st Edition (1996)

Fundamentals of Enzyme Kinetics 4th Edition (2012)

Theory of the Origin, Evolution, and Nature of Life by Andrulis

 Ars Technica has an interesting article that I can't avoid bringing up here. The title of the article is:

"How the craziest f#@!ing "theory of everything" got published and promoted"

The Ars Technica article describes a paper (Theory of the Origin, Evolution, and Nature of Life) published by an assistant professor from Case Western. The author of the paper describes a theory of everything which because of a press release from Case Western manages to get amplified out of all proportion even though the content is highly suspect. Just reading the first sentence is enough to raise a big red flag. That sentence is:

"How life abides by the second law of thermodynamics yet evolutionarily complexifies and maintains its intrinsic order is a fundamental mystery in physics, chemistry, and biology [1]."

There is no mystery here as the author suggests. If he had bothered to read up on Prigogine's and Nichols well-known work on non-equilibrium thermodynamics published decades ago he would have an explanation for this "mystery".

Front-Loading of Flux Control

Originally Posted on December 20, 2019 by hsauro

There is an interesting property that straight chain pathways have which I call front-loading. The phenomenon has been known for some time and I describe it in detail in my recent textbook on Metabolic Control Analysis. Let’s say we have a straight chain of reactions, something like:

  

There are no long-range feedback loops (other than short-range product inhibition) and it is assumed that for a given reaction, increases in substrate cause the reaction rate to increase and increases in the product causes the reaction to decrease. Each reaction is reversible using a simple reversible mass-action rate law: . I assume that  and  are fixed species. Given suitable rate constants and equilibrium connstants for each reaction, the system will admit a steady state. A question to ask is what is the distribution of flux control across the chain? There are two ways to answer this, we can use theory, this gives us what the distribution is and insight into why or we can do simulation but this will just tell us what the distribution is likely to be. For now let’s do a simulation.

The code is shown below. The code defines a four-step straight-chain, runs 10,000 simulations, randomizing the rate constants for each simulation. Each time it does a simulation we compute the flux control at each step, store this information and then at the end, plot histograms of the distribution of the flux control coefficients. Scroll down to see the plots. The average values for the flux control coefficients is shown in the histogram:

The distribution of flux control is shown below. The vertical axis is the frequency and the x-axis the value of the flux control coefficient from zero to 1.0. The red curve corresponds to the flux control for the first step, on average this has the highest control. The yellow distribution corresponds to the last step of the pathway, you can see it is least likely to have any significant flux control. The conclusion is that given a straight chain with a random set of parameters, on average the first step will have the highest control, progressively decreasing as we work our way down the pathway. Why is this? The explanation is a thermodynamic one which in turn boils down to have easy it is for perturbation to travel up and down the chain. So long as the thermodynamic gradient is from left to right, it is easier for perturbations to propagate downstream compared to upstream. Since flux control is really a measure of how much a perturbation has on the steady-state flux, the easier a perturbation can travel the higher the flux control.

This is not to say that it isn’t possible to get excellent flux control in downstream steps, it’s just that fewer parameter combinations will achieve that.

# Monte Carlo simulation of a straight chain pathway
# Samples parameter values while keeping Keq constant
# Plots the distribution of control coefficients
# Not the most elegant code but I wrote it quickly
 
import tellurium as te
import roadrunner
import random
import pylab as pl
 
r = te.loada("""
   J1: $Xo -> S1; k1*(Xo - S1/Keq1);
       S1 -> S2;  k2*(S1 - S2/Keq2);
       S2 -> S3;  k3*(S2 - S3/Keq3);
       S3 -> $X1; k4*(S3 - X1/Keq4);
     
     k1 = 0.1; k2 = 0.1;
     k3 = 0.1; k4 = 0.1;
     
     Keq1 = 4;
     Keq2 = 3;
     Keq3 = 2;
     Keq4 = 1;
     Xo = 5;
     X1 = 0.1;
""")
 
m = r.simulate (0, 10, 100);
 
aC1 = 0; aC2 = 0; aC3 = 0; aC4 = 0;
aC1a = []; aC2a = []; aC3a = []; aC4a = [];
n = 1000
a = 0
upperLimitK = 50
for i in range (0,n):
    r.setValue ('k1', random.uniform(0, upperLimitK))
    r.setValue ('k2', random.uniform(0, upperLimitK))
    r.setValue ('k3', random.uniform(0, upperLimitK))
    r.setValue ('k4', random.uniform(0, upperLimitK))
    
    try:
      r.simulate()    
      r.steadyState()
      C1 = r.getCC ('J1', 'k1')
      C2 = r.getCC ('J1', 'k2')
      C3 = r.getCC ('J1', 'k3')
      C4 = r.getCC ('J1', 'k4')
      aC1 = aC1 + C1
      aC2 = aC2 + C2
      aC3 = aC3 + C3
      aC4 = aC4 + C4
      aC1a.append (C1)
      aC2a.append (C2)
      aC3a.append (C3)
      aC4a.append (C4)
    except:
        a = a + 1
 
print (aC1/n, aC2/n, aC3/n, aC4/n) 
 
bins = 100
 
pl.hist(aC1a, bins=bins, histtype='stepfilled', density=True, color='r', alpha=0.5, label='C1')
pl.hist(aC2a, bins=bins, histtype='stepfilled', density=True, color='b', alpha=0.5, label='C2')
pl.hist(aC3a, bins=bins, histtype='stepfilled', density=True, color='g', alpha=0.5, label='C3')
pl.hist(aC4a, bins=bins, histtype='stepfilled', density=True, color='y', alpha=0.5, label='C4')
print ("fails = ", a)

Stream Plots Using Tellurium

Originally posted on December 18, 2019 by hsauro

Every wanted to do a phase plot for a toggle switch or any two-dimensional model? Here is an example of using a stream plot combined with SBML models using Tellurium based on a simple toggle switch.

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

import numpy as np import matplotlib.pyplot as plt import tellurium as te   Y, X = np.mgrid[0:8:200j, 0:8:200j]   r = te.loada(''' J1: -> x; k1/(1+y^n1) - k2*x; J2: -> y; k3/(1+x^n2) - k4*y;   x = 4; y = -4; k1 = 12; k3 = 12; k2 = 2; k4 = 2 n1 = 4; n2 = 4 ''')   m = r.simulate (0, 8, 100)   U = np.empty_like(X) V = np.empty_like(X)   for (i,j), value in np.ndenumerate(X):     r.x = X[i][j]; r.y = Y[i][j]     U[i][j] = r.J1; V[i][j] = r.J2   plt.figure(figsize=(10,7)) plt.xlabel('x', fontsize='16') plt.ylabel('y', fontsize='16') plt.ylim((-0.1,8)) plt.xlim((-0.1,8)) plt.streamplot(X, Y, U, V, density=[2, 2])   r.steadyStateSolver.allow_presimulation = False   r.x = 5.999; r.y = 0.00462 r.steadyState() plt.plot(r.x, r.y, 'ro', markersize=12)   r.x = 0.0046; r.y = 5.999 r.steadyState() plt.plot(r.x, r.y, 'bo', markersize=12)   r.x = 1.359; r.y = 1.359 r.steadyState() plt.plot(r.x, r.y, 'go', markersize=12)   plt.show()

The colored circles mark the three steady states. The green circle in the middle is on a saddle-node which is unstable (unless you’re exactly on the saddle ridge) and the two on either side are stable nodes.