Tuesday, August 16, 2022

Plot of 1/x using pgfplots

I noticed I couldn't find pgfplots version of the $1/x$ plot that describes $e$. I needed one, so here is the LaTeX I used.



ytick = {1,2},
xtick = {0,1,2,3},
xlabel=$x$, ylabel=$y$,
ymax=2.5, xmin=0, ymin=0,
xlabel style = {anchor=north east},
ylabel style = {anchor=north east}
\addplot [color = blue, name path=A,domain=0:4, line width = 1.4pt, samples=200] {1/x};
\path [name path=B]
    (\pgfkeysvalueof{/pgfplots/xmin},0) --
\addplot [blue!20]
    fill between [of=A and B,
                  soft clip={(1,0) rectangle (2.71,25)},];
\node at (axis cs:1.5,1.2) {$y=1/x$};
\node at (axis cs:1.6,0.31) {Area = $1$};
\node at (axis cs:2.71, 0.82) {$e$};
\addplot[-latex] coordinates
           {(2.71,0.75) (2.71,0.45)};


Friday, July 29, 2022

To paraphrase an old saying....

 You can lead a modeler to a technical solution, but you can't make them use it.

This rephrasing of an old saying was inspired by our reproducibility work. There exist technical solutions to publishing reproducible results, and yet most studies are still not reproducible because we choose not to use those solutions. 

Sunday, June 26, 2022

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:

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. 

Tuesday, April 19, 2022

Lorenz Attractor

I wanted to test out Google's skia 2D library so I wrote a simple interactive Lorenz attractor app over the weekend. Source code and binaries at 


Binaries are only for Windows at the moment but hope to have a Mac version in a couple of weeks.

I used Object Pascal to write it but it should be easily translatable to something like C# and WinForms which were modeled after Object Pascal and the VCL. Looks like skia does a better job at antialiasing lines and also it's quite fast.

I use a simple Euler integration scheme to solve the Lorenz ODEs, nothing fancy but it seems to work perfectly well. 

Tuesday, April 5, 2022

Analogmachine reborn

I decided to rebuild my blog on blogger. Originally I used my own WordPress site but discovered that maintaining it was quite time-consuming. The main problem was stopping the site from being hacked and corrupted.  This happened again recently and it prompted me to move to a more resilient and less demanding platform.

Monday, April 4, 2022

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".

Testing Code Formating using hightlight.js

# Python Program to find the area of triangle

if x == True:
a = 5
b = 6
c = 7

# Uncomment below to take inputs from the user
# a = float(input('Enter first side: '))
# b = float(input('Enter second side: '))
# c = float(input('Enter third side: '))

# calculate the semi-perimeter
s = (a + b + c) / 2

# calculate the area
area = (s*(s-a)*(s-b)*(s-c)) ** 0.5
print('The area of the triangle is %0.2f' %area)

Sunday, April 4, 2021

A note on ex nihilo

In terms of the origin of the universe, "Out of nothing" is something we often hear from theists who argue that something, i.e the universe,  cannot come from nothing from which they conclude that it must have come from something, which they call God. 

The notion of nothing 'before' the universe however is problematic because what is 'nothing'? It is better to refer to the 'nothing' as something that isn't space-time, rather than a simple negation of space-time. What 'out of nothing' seems to imply is that space-time cannot come from space-time which seems reasonable.

All our notations of cause and effect and the logical absolutes only operate within our own space-time. It is impossible for us to contemplate anything that isn't space-time, so asking questions about creation is extremely difficult if not impossible. Jumping to God as an explanation, however, makes no sense since we can't comprehend anything outside space-time, let alone an all-powerful sentient being to which the bible attaches all sorts of very specific properties. 

Thursday, September 3, 2020

Mathpix Snip

Originally Posted on  by hsauro

I came across this amazing tool that can convert images of math equations into LaTeX format. The tool can be found at https://mathpix.com/. I’ve tried it on a number of texts including some not so clear and it does a fantastic job of converting to LaTeX. Here is a screen showing part of a page from Paul’s Online Notes:

The way it works is you select the screen icon on the mathpix tool, the entire screen goes black and white, then we draw a square around the section we want to convert and that’s it. In this case, it generates the following latex

We’ll start with finding the derivative of the sine function. To do this we will need to use the definition of the derivative. It’s been a while since we’ve had to use this, but sometimes there just isn’t anything we can do about it. Here is the definition of the derivative for the sine function.
\frac{d}{d x}(\sin (x))=\lim _{h \rightarrow 0} \frac{\sin (x+h)-\sin (x)}{h}
since we can’t just plug in $h=0$ to evaluate the limit we will need to use the following trig formula on the first sine in the numerator.
\sin (x+h)=\sin (x) \cos (h)+\cos (x) \sin (h)
Doing this gives us,
\frac{d}{d x}(\sin (x)) &=\lim _{h \rightarrow 0} \frac{\sin (x) \cos (h)+\cos (x) \sin (h)-\sin (x)}{h} \\
&=\lim _{h \rightarrow 0} \frac{\sin (x)(\cos (h)-1)+\cos (x) \sin (h)}{h} \\
&=\lim _{h \rightarrow 0} \sin (x) \frac{\cos (h)-1}{h}+\lim _{h \rightarrow 0} \cos (x) \frac{\sin (h)}{h}
As you can see upon using the trig formula we can combine the first and third term and then factor a sine out of that. We can then break up the fraction into two pieces, both of which can be dealt with separately.

Which when processed by LaTeX becomes:

This is rendered inside the mathpix tool but you’ll notice there isn’t a significant difference between the original and the converted image. I’ve converted some fairly rough images and it generally succeeds. It also gives you a confidence level on how well it thinks it’s done. It took under a second to generate the LaTeX.

As a harder test, I decided to attempt to translate a page from Jim Burns’ thesis. This is a thesis from the 1970s that was typed and the equations a combination of typed characters and hand drawn. The following image shows page 93 which is part of the proof for the connectivity theorem.

And here is the image analyzed by mathpix. Remarkably the conversion is almost perfect, the equations, in particular, are translated almost without error, even getting the subscripts on the subscripts correct. It got the delta F1 wrong at the start and it interpreted a mark on the paper as an apostrophe. I tried other pages that included derivatives and these converted without incident.

Wednesday, April 22, 2020

PID Control Demonstration

Those of you who are interested in feedback control may have come across the infamous PID controller. To many engineers, PID control is the answer to all their problems. What is it for? PID controllers are a sophisticated way to stabilize a system that is under continual threat of disturbance. The simplest example is the cruise control in a car. Most of us probably set the cruise control and forget about it. What’s happening behind the scenes is however much more interesting. As we go up a hill we notice the engine revving up to maintain the set cruise speed and when we go down a hill the engine slows down. The simplest feedback mechanism is where the controller compares the speed of the car with the set cruise speed and adjusts the engine accordingly. The secret to a good controller such as this is what it does with the error it finds between the actual speed of the car and the desired speed. The simplest solution is to adjust the speed of the engine directly in proportion to the difference, this is called proportional control and is most familiar to engineers. The problem with proportional control is twofold:

1. Too much proportional control and the response overshoots, perhaps so much so that it goes into uncontrolled oscillations, not a good thing in a car.

2. You usually have to set the cruise control slightly higher than the speed you actually want because a proportional controller finds it difficult to actually match the set speed unless the proportional control is strong enough (See problem 1!)

Both these problems can be solved by combining proportional control (The P in PID) with two others, integral (The I in PID) and derivative control (The D in PID). Derivative control is perhaps the easiest to understand. What derivative control does is measure how fast the error between desired and actual engine speed changes. If the difference is changing too fast then we want to slow down the response to avoid overshooting. Derivative control will therefore help get rid of the instability (Too much derivative control can however also result in instability).

The second problem is that the speed of the car doesn’t actually settle to the desired speed. This is fixed by using integral control. This is where the error is continually added up to form a net error and it is the net error than is then used to control the speed of the engine. So long as there is even a little bit of error (ie Desired speed 50 mph, actual speed 49 mph), the net error will increase and drive the car speed to the desired speed until the actual error is zero. Integral is particularly useful if there is a continual disturbance at work, for example, a headwind.

To demonstrate these capabilities, I have written a small PID control demonstration tool (for windows only I am afraid) that simulates each of the three elements in a PID controller. The screenshot below shows the application. A user can control the strength of the different controlling elements or remove them all together. To test the response a user can apply a step change to the desired car speed and see how the system responds. The simulation is continuous so it is easy to see the effects.

Partly supported by the NSF

Friday, December 20, 2019

Front-Loading of Flux Control

Originally Posted on  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:

  \[ X_o \rightleftharpoons S_1  \rightleftharpoons S_2 \rightleftharpoons \hdots \rightleftharpoons  S_n \rightleftharpoons X_1 \]

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: k_1 (S - P/K_{eq}). I assume that X_o and X_1 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 few 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))
      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)
        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)