top of page

The Stoichiometry Matrix in Reaction and Epidemiological Models

Networks of chemical or biological reactions can be represented by systems of coupled differential equations. Solving these provides the concentrations of all species in time. Compartmental epidemiological models can be represented in the same way, with the solutions providing the number of people in each epidemiological state as a function of time. However, it is tedious to set up and solve a set of coupled differential equations every time we want to investigate a reaction or epidemiological model, so a general tool would take a lot of work out. This is provided by the stoichiometry matrix. By being able to cast a chemical, biological or epidemiological problem in terms of its stoichiometry matrix, we can use a single formalism (and so, single computer program and data fitting system) for any network of arbitrary complexity.

To see how we get to a general expression for a reaction network or epidemiological model, start with an example:

In this reaction network, A and B react to transform A into another B irreversibly. B then has two transformation pathways into either C or D (One single B transforms to one single C or D). So B is a catalyst for its own synthesis, but is also subject to transformation. If you are more interested in epidemiology than chemistry or systems biology, replace "A,B,C,D" with "S,I,R,D" as this network represents the "Susceptible-Infectious-Recovered" model including a mortality rate with D representing Deceased state. The system of reaction equations can be written as:

A,B,C,D denote time-dependent species concentrations or counts. The dot represents multiplication. Now, re-write the reaction equation for the consumption of some species as:

Here, S is the number of C_i molecules involved (its stoichiometry) in the particular reaction, k the rate constant for the reaction, and P_{m} the reaction orders for each species C_m involved. So the rate at which C_i is taken out of the system is its stoichiometry, multiplied by its rate constant, multiplied by the product of all species involved raised to the powers of their reaction orders (0 for species not involved). However, there will be many C_i, and many reactions involved. We can consider each "point" in our set of reaction equations to be a node in a network. In our example, there are four nodes. A rate constant specifies communication rates between nodes. Forward and reverse traffic must be separately considered if present (it is not in our first example), with every individual reaction arrow therefore a communication channel with and associated communication rate given by its rate constant. The generalisation of the above equation for all species follows as:

We have now introduced a matrix of stoichiometries S. Each row indexes a species, each column indexes a reaction (or network channel). The element S_{ij} is the number of molecules of species i involved in reaction j. The element is negative if molecules of species indexed i are taken out of the system by channel j (a "forward reaction"), positive if they are put into the system (a "reverse reaction"). The reaction order matrix is similarly arranged, but with different elements. Its elements are the reaction orders for each species, taking the value zero if species m is not involved. The element P_{mj} is the reaction order for species m and channel j. The reaction order matrix ensures that only the relevant species are retained in the product. The stoichiometry matrix and reaction order matrix together connect the state space to the network channel space. The rate constants are arranged as a vector with k_j containing the rate constant for channel j. Note that we may express this in vector-matrix form as:


Now, how might we set up S and P for our network? Channels and species have first to be enumerated. The ordering is our choice, but a sensible ordering here is:

for the ordering of time-dependent species concentrations (or counts) and

for the vector of rate constants. This sets the ordering of the species (states c(t)) and channels (rate constants k). To assemble row one the stoichiometry matrix (for species A), we identify the number of molecules of A associated with each of the channels in k in vector form to give:

So the k1 channel takes 1 A out of the system, and neither of the other channels communicate with A. Following that lead, it is easy enough to complete the stoichiometry matrix:

Recall from the definitions of c(t) and k that the rows are ordered A,B,C,D (or S,I,R,D if you want an epidemiological interpretation), and the columns k1,k2,k3. We now need the reaction orders matrix P. In this system, the kinetics are either first-order if a species is involved or zero-order if not. So P may be constructed simply as:

So looking at reaction 1, represented by column 1, it is 1st order in A and B, and so on.


A second example: An enzyme-catalysed reaction

As another simple example of setting up a system of rate equations for a reaction network, take the most basic model of a enzyme-catalysed reaction. Here, enzyme E and substrate S reversibly form complex ES, which reversibly forms E and product P, with all reactions first order with respect to all reagents. The reaction equations are:

We choose to order the channels:

And we choose to order the species:

We will start with the first row of the stoichiometry matrix:

So the k_{12} channel takes one E out of the system, k_{21} puts one E into the system, k_{23} puts one E into the system and k_{32} takes one E out of the system. Following that lead, it is easy enough to complete the stoichiometry matrix:

In this system, the kinetics are either first-order if a species is involved or zero-order if not. So P may be constructed simply as:

It should be easy enough to multiply these out and verify the rate equations for this system as:



bottom of page