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

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 

https://pypi.org/project/sphinx-execute-code-python3/ 

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.

A simple linear pathway

Analysis of a simple pathway:

This model is model a) from the paper by Tyson et al (2003) "Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell". I changed their R to P however so I could use the symbol R for response.

Assume $ v_1 = k_1 S $ and $ v_2 = k_2 P $ Then $\frac{dP}{dt} = k_1 S - k_2 P $. At steady-state this equals 0 so that $$ P = \frac{k_1 S}{k_2} $$ In other words P is a linear function of S. In terms of logarithmic senstivities or the reponse coefficient: $$ R^P_S = \frac{dP}{dS} \frac{S}{P} = \frac{k_1}{k_2} \frac{k_2}{k_1} = 1 $$ That is a 1% change in S will lead to a 1% change in the steady-state level of P. This is important because in their pathway (d, shown below) they use this property to ensure that when S also regulates a secondary pathway that in turn up regulates the consumption step v2, the level of P remains unchanged.

This property is not very robust however and any change to the rate law on $v_1$ will result in the pathway failing to maintain P constant. For example if $v_1 = k_o + k_1 S$, where we've added a basal rate of $k_o$, then the response becomes: $$ R^P_S = = \frac{k_1 S}{k_o + k_1 S} $$ that is the response is no longer proportional. Change S in this case will result in P changing. The code below will model this system. One important point to make is that this isn't an example of integral control.

import tellurium as te
import matplotlib.pyplot as plt
  
r1 = te.loada("""
     -> P; k0 + k1*S
     P ->; k2*X*P
     -> X; k3*S
     X ->; k4*X
     
     k1 = 1; k2 = 1
     k3 = 1; k4 = 1
     
     # Set basal rate to zero 
     k0 = 0; S = 1
     
     # Change signal, P won't change
     at time > 10: S = S*2
     
     # Change basal rate and set S back to what it was
     at time > 25: k0 = 0.3, S = 1;
     
     # Change signal, this time P will change
     at time > 40: S = S*2      
""")

m = r1.simulate(0, 60, 200, ['time', 'P', 'S'])
plt.plot (m['time'], m['P'], label='P')
plt.plot (m['time'], m['S'], label='S')
plt.text(2, 0.75, "Basal = 0")
plt.text(9, 0.9, "Change S")
plt.text(14, 1.2, "P restored")
plt.text(20, 0.75, "Set basal > 0")
plt.text(35, 0.9, "Change S, P not restored")
plt.legend()

Plotting a steady-state plot

It's always surprising to find students who don't quite get what a steady-state is. Here is a simulation of a three step pathway that might help. You'll need to install tellurium for this to work.

$X_o \stackrel{v_1}{\rightarrow} S_1 \stackrel{v_2}{\rightarrow} S_2 \stackrel{v_3}{\rightarrow}$

we assume that $X_o$ is fixed and does not change in time. I don't care where $S_3$ goes.

If we start with zero concentrations for $S_1$ and $S_2$ we get the following plots. The left plot shows the change in concentration and the right plot the reaction rates. Note that all three reaction rates approach the same rate since at steady-state all rates must be equal.

Here is the python code that will generate that plot. To get the reaction rate symbols above the reaction arrows I used a bit of a hack by printing two text strings, one above the other. For some reason I couldn't get matplotlib to recognise the LaTeX command stackrel.

import tellurium as te
import matplotlib.pyplot as plt

r = te.loada("""
    J1: $Xo -> S1; k1*Xo - k2*S1
    J2:  S1 -> S2; k3*S1 - k4*S2
    J3:  S2 -> ;   k5*S2
      
      k1 = 0.1; k2 = 0.04
      k3 = 0.14; k4 = 0.09
      k5 = 0.16
      Xo = 10
""")

m = r.simulate (0, 60, 100)
plt.subplot (1,2, 1)
plt.plot(m['time'], m['[S1]'], label='S1')
plt.plot(m['time'], m['[S2]'], label='S2')
plt.ylim((0, 9))
plt.xlabel('Time')
plt.ylabel('Concentration')
plt.text(5, 8.4, 'Approach to steady-state')
plt.text(5, 7.8, '$X_o$ is fixed')
plt.text(14, 1.5, r'$\quad\  v_1 \quad v_2 \quad v_3$', fontsize = 12) 
plt.text(14, 1, r'$X_o \rightarrow S_1 \rightarrow S_2 \rightarrow$', fontsize = 12) 
plt.legend(bbox_to_anchor=(0.4, 0.5))

# Next generate the reaction rates
r.reset()
m = r.simulate (0, 40, 100, ['time', 'J1', 'J2', 'J3'])
plt.subplot (1,2, 2)
plt.plot(m['time'], m['J1'], label='$v_1$')
plt.plot(m['time'], m['J2'], label='$v_2$')
plt.plot(m['time'], m['J3'], label='$v_3$')
plt.xlabel('Time')
plt.ylabel('Flux')
plt.text(14, 0.2, r'Reaction rates', fontsize = 12) 
plt.legend()
plt.tight_layout()
plt.show()

1888 Math Exam

Nothing to do with cells or modeling, but I recently purchased a 1885 copy of John Casey's rendition of Euclid's Elements. John Casey was born in Limerick, Ireland, and became a lecturer in mathematics at University College Dublin. He wrote one of the more well-known editions of Euclid's elements. In his preface he states:

"This edition of the Elements of Euclid, undertaken at the request of the principals of some of the leading Colleges and Schools of Ireland, is intended to supply a want much felt by teachers at the present day—the production of a work which, while giving the unrivalled [sic] original in all its integrity, would also contain the modern conceptions and developments of the portion of Geometry over which the Elements extend."

You can find a copy at Project Gutenberg (https://www.gutenberg.org/ebooks/21076) if you are curious. 

What's interesting about the book I received, is not so much the geometry, but what I found inside. 

Inside was a copy of a math Pass Examination paper from 1888, presumably from University College Dublin but at least somewhere in Ireland (I purchased the book from Dublin). Update: There are two inscriptions in the book. The first is a simple "W.H Dunlop May 1886".  The second is more interesting and written in a very elegant cursive style where the book was transferred two years later to: Michael J Buckley. Catholic University Dublin. 19.1.'88. This is the name of the first university dedicated to accepting Catholics from Ireland (Trinity College was the Anglican University founded by Elizabeth I, so you can imagine the problem). However, the university didn't do too well financially and had problems awarding degrees because it didn't have a royal charter. In 1908/1909, it became the University College Dublin with its own charter. 

The exam had 10 questions, one involving reaping a field, another about selling a horse, and another about carpeting a room. The remaining 7 are pure math questions, with a number of them being numerical estimation questions, eg, roots, squares etc. I think modern students could do these questions, though they might struggle with the numeric computations unless they have a calculator at hand. I took a photograph of the exam for all to see:

Optimal distribution of enzymes that maximizes flux

Here is an interesting metabolic question to ask. Consider a pathway with two steps: $$\large X_0 \stackrel{e_1}{\rightarrow} S \stackrel{e_2}{\rightarrow} $$ where the first step is catalyzed by enzyme $e_1$ and the second step by $e_2$. We can assume that a bacterial cell has a finite volume, meaning there is only so much space to fit all the proteins necessary for the bacterium to live. If the cell make more of one protein, it presumably needs to find the extra space by down expressing other proteins. If the two step pathway is our cell, this translates to there being a fixed amount of protein that can be distributed betwen the two, that is: $$ e_1 + e_2 = e_t $$ Presumably evolutionary pressure will maximize the flux through the pathway per unit protein. This means that given a fixed total amount of protein $e_t$ there must be a particular distribution of protein between $e_1$ and $e_2$ that maximizes flux. For example, if most of the protein were allocated to $e_1$ and very little to $e_2$ then the flux would be quite low. Likewise if most of the protein were allocated to $e_2$ and very little to $e_1$ then again the flux would be low. Between these two extremes there must be a maximum flux. To show this is the case we can do a simple simulation shown in the python code below. This code varies the levels of $e_1$ and $e_2$ in a loop such that their total is always 1 unit. Each time we change $e_1$ we must change $e_2$ by the same amount in the opposite direction in order to keep the constant fixed. As we do this we collect the steady-state flux and the current level of $e_1$. Finally we plot the flux versus $e_1$.

import tellurium as te
import roadrunner
import matplotlib.pyplot as plt

r = te.loada("""
     J1: $Xo -> S; e1*(k1*Xo - k2*S)
     J2: S ->; e2*k3*S
     
     Xo = 1
     k1 = 0.5; k2 = 0.2
     k3 = 0.45
     e1 = 0.01; e2 = 0.99
""")

r.steadyState()

x = []; y = []; 
for i in range (49):
    r.e1 = r.e1 + 0.02
    r.e2 = 1 - r.e1  # Total assumed to be one
    r.simulate()
    r.steadyState()
    x.append (r.e1)
    y.append (r.J1)
    
plt.figure(figsize=(12,10))
plt.grid(b=True, which='major', color='#666666', linestyle='-')
plt.grid(b=True, which='minor', color='#999999', linestyle='-', alpha=0.2)
plt.minorticks_on()
plt.plot (x, y, linewidth = 2)
plt.xlabel('$e_1$', fontsize=36)
plt.ylabel('Flux', fontsize=36)

 

The following graph show the results from the simulation (the font sizes might be on the big size if you have a low res monitor):

You can see the flux reaching a maximum at around $e_1 = 0.6$, meaning also that $e_2 = 0.4$ in order to keep the total fixed. The actual position of the peak wil depend on the rate constants in the rate expressions. Let's assume we are at the maxium. The slope, $dJ/de_1$, at the maxium is obviously zero. That means if we were to move a small amount of protein from $e_1$ to $e_2$ the flux won't change. We can write this experiment in terms of the two flux control coefficients: $$ \frac{\delta J}{J} = 0 = C^J_{e_1} \frac{\delta e_1}{e_1} + C^J_{e_2} \frac{\delta e_2}{e_2} $$ However, we know that in this particlar experiment the change in $e_1$ is the same but oppisite to the change in $e_2$. That is $\delta e_1 + \delta e_2 = 0$ or $$ \delta e_1 = -\delta e_2$$ Replacing $e_2$ with $-\delta e_1$ gives: $$ \frac{\delta J}{J} = 0 = C^J_{e_1} \frac{\delta e_1}{e_1} - C^J_{e_2} \frac{\delta e_1}{e_2} $$ $$ \frac{\delta J}{J} = 0 = C^J_{e_1} \frac{1}{e_1} - C^J_{e_2} \frac{1}{e_2} $$ $$ C^J_{e_1} \frac{1}{e_1} = C^J_{e_2} \frac{1}{e_2} $$ Giving the final result: $$ \frac{C^J_{e_1}}{C^J_{e_2}} = \frac{e_1}{e_2} $$ That is, when the protein distribution is optimized to maximize the flux, the flux control coefficients are in the same ratio as the ratio of enzyme amounts. This generalizes to any size pathway with multiple enzymes. This gives us a tantgilizing suggestion that we can obtain the flux control coeficients just by measuring the protein levels.

There is obviously a lot more one can write there and maybe I do that in future blogs but for now you can get further information:

S. Waley, “A note on the kinetics of multi-enzyme systems,” Biochemical Journal, vol. 91, no. 3, p. 514, 1964.

J Burns: “Studies on complex enzyme system.” https://era.ed.ac.uk/handle/1842/13276, 1971 (page 141-) I have a LaTeX version at: https://github.com/hsauro/JumBurnsThesis

Guy Brown, Total cell protein concentration as an evolutionary constraint on the metabolic control distribution in cells,” Journal of theoretical biology, vol. 153, no. 2, pp. 195–203, 1991.

E. Klipp and R. Heinrich, “Competition for enzymes in metabolic pathways:: Implications for optimal distributions of enzyme concentrations and for the distribution of flux control,” Biosystems, vol. 54, no. 1-2, pp. 1–14, 1999

Sauro HM Systems Biology: An Introduction to Metabolic Control Analysis, 2018

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.

How not to Comment Code

I recently came across one of the ten simple rules aricles in PLoS Comp Bio:

"Ten simple rules for tackling your first mathematical models: A guide for graduate students by graduate students" by Korryn Bodner et al

https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1008539

One thing that struck me was Rule 5 on coding best practices with commenting being one of the discussion points. What struck me was their screen shot of a documented function shown below (in R):

My take on commenting is that it should be used to add human readable metadata on elements of a program that are not immediately obvious.

Most of the time, code should be sufficiently readable to indicate what it's doing. Obviously some languages are better than others when desribing an algorithm but it is also dependent on the programmer. I've seen code written in clear languages that are unintelliglbe, but I've also seen code written in poorly expressible languages that are easily readable. Although the programming language itself can influence code reability I think the programmer has much more influence.

But back to Rule 5. In the example you'll see something like:

# calculate the mean of the data

u <- mean (x)

This is completely redundant, as the coding states what it is going to do. In fact the authors comment every line like this. If anythng, I think the extent of comments actually hinders the reabilty of the code. The code itself is mostly clear as to what it is doing. There may be a justification to include a comment on next line that computes the standard error because the variables names are so badly chosen, e.g what does the following line do:

s <- sd(x)

sd might stand for standard deviation but the rest of the line offers no clue. If it had been written as:

standardDeviation <- sd(x)

It would have been much clearer, instead the authors add a comment to make up for poor choice of variable names. They also give the function itself a nondescriptive name, in this case ci. It would have been better to write the function using getConfidenceInterval or similar:

getConfidenceInterval <- f (data) {

etc

Euclid's Elements

Nothing to do with cells and networks, but I have an interest in Euclid's Elements and decided to publish a new edition of Book I. The difference with previous editions is that this one is in color and also has a chapter on the history of the elements, together with commentaries on each proposition. For those unfamiliar with Euclid's Elements, it's a series of books (more like chapters in today's language) laying the foundation for geometry and number theory. Book I focuses on the foundation of geometry culminating in a proof for Pythagoras' theorem and other important but less known results. The key innovation is that it describes a deductive approach to geometry. It starts with definitions and axioms from which all results are derived using logical proofs.

It occurred to me that something similar could be done with deriving the properties of biochemical networks. For example, we might define the following three primitives:

I. Species

II. Reaction

III. Steady-state

We might then define the following axioms:

I. A species has associated with it a value called the concentration, x_i.

II. All concentrations are positive.

III. A reaction has a value associated with it called the reaction rate, v_i.

IV. Reaction rates can be negative, zero, or positive.

V. A reaction transforms one or more species (reactants) into one or more other species (products).

VI. The reaction rate is a continuous function of the reactants and products.

VII. The rate of change of a species can be described using a differential equation, dx/dt

VIII. All steps are reversible unless otherwise stated (may this can be derived?)

etc

Given these axioms, we could build a series of propositions. This might be an interesting exercise to do. Some of the more obvious propositions would be the results from metabolic control analysis, such as the summation and connectivity theorems.

Link to Amazon

Branched Pathway in TikZ

 I couldn't resist, here are some simple branched pathways using tikz. Note you can change the vertical splay and horizontal distance for the nodes by changing the node distance in the argument to tikzpicture. I also added some colors wshich really don't match, someone with more artistic tallent could do better.

\documentclass[12pt]{article}

\usepackage{tikz}
\usetikzlibrary{positioning}
\usepackage{xcolor}

\begin{document}

\begin{tikzpicture}[node distance=0.4cm and 1cm]
\node (Xo) at (0,0) {};
\node[right=of Xo] (S1) {\small $S_1$};
\node[above right=of S1] (S2) {};
\node[below right=of S1] (S3) {};

\draw[-latex, thick, blue] (Xo) -- (S1) node[pos=0.5,above,gray] {$v_1$};

\draw[-latex, thick, orange] (S1) -- (S2) node[above=0.15, left,pos=0.5,gray] {$v_2$};
\draw[-latex, thick, green!45!black] (S1) -- (S3) node[below=0.15, left,pos=0.5,gray] {$v_3$};
\end{tikzpicture}

\vspace{1cm}

\begin{tikzpicture}[node distance=0.4cm and 1cm]
%\draw [help lines,step=.1] (0,-6) grid (6,6);
\node (Xo) at (0,0) {};
\node[right=of Xo,yellow!20!red] (S1) {\small $S_1$};
\node[above right=of S1,blue!80!red] (S2) {\small $S_2$};
\node[below right=of S1,green!80!red] (S3) {\small $S_3$};

\node[above right=of S3] (S4) {};
\node[below right=of S3] (S5) {};

\draw[-latex, thick, blue] (Xo) -- (S1) node[pos=0.5,above,purple] {$v_1$};
\draw[-latex, thick, blue] (S1) -- (S2) node[above=0.15, left,pos=0.5,purple] {$v_2$};
\draw[-latex, thick, blue] (S1) -- (S3) node[below=0.15, left,pos=0.5,purple] {$v_3$};

\draw[-latex, thick, blue] (S3) -- (S4) node[below=-0.12, left,pos=0.5,purple] {$v_4$};
\draw[-latex, thick, blue] (S3) -- (S5) node[below=0.05, left,pos=0.5,purple] {$v_5$};
\end{tikzpicture}

\end{document}

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