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.

Many years ago I wrote a windows application that could be used to build biochemical models visually. Origially called JDesigner but subsequently renamed to PathwayDesigner. Over the years the code grew larger and somewhat unwieldy. In addition it couldn't run on the Mac. I have decided therefore to rewrite it using more modern tech so that is can run on the Mac as well as Windows. Initially it won't include simulation capabilities since simulators are dime a dozen. It also won't support SBGN or the SBML layout/render standard since htey are hard to support. The GitHub repo can be found here, currenlty at version 0.1.1 https://github.com/hsauro/PathwayDesigner

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:

A small editor App to automatically format JSON files

I've been very remiss in my postings. Grant proposals and administration consume a huge amount of time. I estimate paperwork now takes up between 4 and 5 months a year of my time. ChatGPT is starting to help with some mundane tasks and hopefully, that will increase over time. Ideally, I’d like to live in a world where all bureaucratic chores are relegated to an AI assistant so we AI humans can focus on more creative activities.

I’m not sure how much you’ve had to work with JSON files, they pop up everywhere and I’ve also ended up using them to store information. The problem I always face is being able to edit them. JSON in its purest form is just a single line of text so that when loaded into your favorite editor you just see one long line of text which isn’t great for editing. Usually, I will paste the text into an online formatter and paste the pretty version of the JSON back into my editor where I can work with it. Yes, I could load the JSON into a specialist editor, but they are even more tedious to use. Yes, I could install a plugin for Notepad++ and try to remember keystrokes to format the JSON or I could figure out how to install yet another plugin in the sublime text and remember yet another keystroke.

I therefore spent probably an hour creating a simple editor that when I load a JSON file it automatically formats the text. I uploaded it to GitHub. More important here is the link to the binaries. Appologies to Mac owners, it ony works on Windows. I don't have a good open source editor component for the Mac platform.

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:

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:

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; 

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