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

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

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 

https://github.com/hsauro/Lorenz_fmx

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. 

Marsaglia’s MWC (multiply with carry) Random Number Generator for Object Pascal

Originally Posted on December 14, 2019 by hsauro 

In the past, I’ve tended to use the Mersenne Twister for generating pseudo-random numbers. This is a widely used generator because it has a very long period of 2^{19937} – 1. The one issue is that it’s not the fastest generator and this is particularly an issue when drawing random numbers for implementing the Gillespie direct method. A much faster variant is the Multiply-with-carry (MWC) random number generator. This uses much simpler arithmetic which means it’s faster. The algorithm can be described in three lines:

  

  

  

The symbol >> means shift the bits to the right by the given number and << is similar but moves the bits to the left. In object pascal, these operations are represented by shr and shl respectively. The generator needs two values to be initialized, m_z and m_w I've recently implemented the Gillespie method in three languages, C, Go and Object Pascal. Object Pascal has good support for the Mersenne Twister but I had to create my own MWC generator. The code below shows the object pascal implementation:

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

// ------------------------------------------------------------ // George Marsaglia's MWC (multiply with carry) generator. // This code translated from C# to Object Pascal and obtained from // https://www.codeproject.com/Articles/25172/Simple-Random-Number-Generation   unit uRandMWC;   interface   Uses System.Types;   // Call this before doing anything procedure initMWC;   // If you want to change the seed call this procedure setSeed (seed : cardinal);   // Call this to get a random number between 0 and 2^32 function  getUint() : cardinal;   // Call this to get a uniformly distributed random number between 0 and 1.0 function getUdouble : double;   implementation   // In Delphi 32 and 64 bit, cardinal is an unsigned 32-bit integer   var   m_z, m_w : cardinal;   procedure setSeed (seed : cardinal); begin   m_z := seed;   m_w := 678934; end;   procedure initMWC; begin   m_z := 467567;   m_w := 125681; end;   function getUint() : cardinal; begin   m_z := 36969 * (m_z and 65535) + (m_z shr 16);   m_w := 18000 * (m_w and 65535) + (m_w shr 16);   result := (m_z shl 16) + m_w; end;   function getUdouble : double; var u : dword; begin   u := getUint();   // The magic number below is 1/(2^32 + 2).   // The result is strictly between 0 and 1.   result := double (u + 1.0) * 2.328306435454494e-10; end;   end.

What about some tests? There are random generator test systems such as TestU01 and dieharder but they are written to run on Linux and I didn't have time to set something up on windows plus the usage docs are not that great. Instead I did three simple sanity tests. The simplest is to check that the mean value for a sample of random numbers is 0.5. Another test is to make sure that, say out of 10,000 random numbers, 10% are between 0 and 0.1, etc. Ideally there should be statistical tests done to check that these values are within certain tolerances. The third test is a visual test where I generate 32768 random numbers and map them on to a bitmap. Each random number is 32 bits long and each bit represents a pixel in the bitmap. 32768 numbers means I generate a bitmap 1024 by 1024 pixels. Thus 32 numbers occupy each row of the bitmap. I wrote a simple Delphi console application that calls the random number generator and creates the bitmap. The main part of this code is shown below. For a comparison I also implemented the NOT to be used randu generator and also the Delphi's built in generator. The approach I take isn't very sophisticated and I am sure it could be made shorter, but for each random number I generate, I convert it into a string of 1's and 0's. I then work my way through the string setting pixels in a bitmap. I do this 32 times for each row in the bitmap.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

program RandomGen;   {$APPTYPE CONSOLE} {$R *.res}   uses   System.SysUtils,   System.Types,   Vcl.Graphics,   uRandMWC in 'uRandMWC.pas';   var r1 : double;       ri : cardinal;     i, j, k : integer; bmp : TBitmap;     sb : string;       iCol : integer;     state : cardinal;  count : integer;   // Get a single bit as position Bit from cardinal aValue function getBit(const aValue: cardinal; const Bit: Byte): Boolean; begin   Result := (aValue and (1 shl Bit)) <> 0; end;   // Converts a cardinal (32-bits) into a string of ones and zeros function getBinaryString (aValue : cardinal) : string; var i : integer; begin   result := '';   for i := 0 to 32 - 1 do       if getBit (aValue, i) then          result := result + '1'       else         result := result + '0' end;   function randu : cardinal; begin    state := (65539 * state) mod 2147483648;    result := state; end;   begin   state := 25; // starting value for radu   writeln ('wait.....');   initMWC();     bmp := TBitmap.Create;   bmp.PixelFormat := pf1bit;   bmp.Width := 1024;   bmp.Height := 1024;   count := 0;   for i := 0 to bmp.Height - 1 do       begin       iCol := 0;       for j := 0 to bmp.Width div 32 - 1 do           begin           // Choose a random number generator           //ri := randu;           ri := getUint();           //ri := random (4294967296-1);           sb := getBinaryString (ri);           for k := 0 to 32 - 1 do               begin               if sb[k+1] = '0' then                  begin                  bmp.Canvas.Pixels[iCol,i] := clWhite;                  end               else                  begin                  bmp.Canvas.Pixels[iCol,i] := clBlack;                  end;               inc (iCol);               end;           end;       inc (count, bmp.Width div 32);       end;     bmp.SaveToFile('image.bmp');   writeln ('count = ', count);   writeln ('Press any key to coninue');   readln; end.

The results of running this code is below. First randu to show you how bad it is, note the patterns and lines in the bitmap:

The next figure shows the results from Delphi’s builtin random number generator, which isn’t so great either:

Finally, the MWC algorithm, which you can see from the image the distribution of pixels appears to show no pattern: