Archive for the ‘dynamical systems’ Category

Models for photosynthesis, part 4

September 19, 2016

In previous posts in this series I introduced some models for the Calvin cycle of photosynthesis and mentioned some open questions concerning them. I have now written a paper where on the one hand I survey a number of mathematical models for the Calvin cycle and the relations between them and on the other hand I am able to provide answers to some of the open questions. One question was that of the definition of the Pettersson model. As I indicated previously this was not clear from the literature. My answer to the question is that this system should be treated as a system of DAE (differential-algebraic equations). In other words it can be written as a set of ODE \dot x=f(x,y) coupled to a set of algebraic equations g(x,y)=0. In general it is not clear that this type of system is locally well-posed. In other words, given a pair (x_0,y_0) with g(x_0,y_0)=0 it is not clear whether there is a solution (x(t),y(t)) of the system, local in time, with x(0)=x_0 and y(0)=y_0. Of course if the partial derivative of g with repect to y is invertible it follows by the implicit function theorem that g(x,y)=0 is locally equivalent to a relation y=h(x) and the original system is equivalent to \dot x=f(x,h(x)). Then local well-posedness is clear. The calculations in the 1988 paper of Pettersson and Ryde-Pettersson indicate that this should be true for the Pettersson model but there are details missing in the paper and I have not (yet) been able to supply these. The conservative strategy is then to stick to the DAE picture. Then we do not have a basis for studying the dynamics but at least we have a well-defined system of equations and it is meaningful to discuss its steady states.

I was able to prove that there are parameter values for which the Pettersson model has at least two distinct positive steady states. In doing this I was helped by an earlier (1987) paper of Pettersson and Ryde-Pettersson. The idea is to shut off the process of storage as starch so as to get a subnetwork. If two steady states can be obtained for this modified system we may be able to get steady states for the original system using the implicit function theorem. There are some more complications but the a key step in the proof is the one just described. So how do we get steady states for the modified system? The idea is to solve many of the equations explicitly so that the problem reduces to a single equation for one unknown, the concentration of DHAP. (When applying the implicit function theorem we have to use a system of two equations for two unknowns.) In the end we are left with a quadratic equation and we can arrange for the coefficients in that equation to have convenient properties by choosing the parameters in the dynamical system suitably. This approach can be put in a wider context using the concept of stoichiometric generators but the proof is not logically dependent on using the theory of those objects.

Having got some information about the Pettersson model we may ask what happens when we go over to the Poolman model. The Poolman model is a system of ODE from the start and so we do not have any conceptual problems in that case. The method of construction of steady states can be adapted rather easily so as to apply to the system of DAE related to the Poolman model (let us call it the reduced Poolman model since it can be expressed as a singular limit of the Poolman model). The result is that there are parameter values for which the reduced Poolman model has at least three steady states. Whether the Poolman model itself can have three steady states is not yet clear since it is not clear whether the transverse eigenvalues (in the sense of GSPT) are all non-zero.

By analogy with known facts the following intuitive picture can be developed. Note, however, that this intuition has not yet been confirmed by proofs. In the picture one of the positive steady states of the Pettersson model is stable and the other unstable. Steady states on the boundary where some concentrations are zero are stable. Under the perturbation from the Pettersson model to the reduced Poolman model an additional stable positive steady state bifurcates from the boundary and joins the other two. This picture may be an oversimplification but I hope that it contains some grain of truth.

NFκB

May 1, 2016

NF\kappaB is a transcription factor, i.e. a protein which can bind to DNA and cause a particular gene to be read more or less often. This means that more or less of a certain protein is produced and this changes the behaviour of the cell. The full name of this transcription factor is ‘nuclear factor, \kappa-light chain enhancer of B cells’. The term ‘nuclear factor’ is clear. The substance is a transcription factor and to bind to DNA it has to enter the nucleus. NF\kappaB is found in a wide variety of different cells and its association with B cells is purely historical. It was found in the lab of David Baltimore during studies of the way in which B cells are activated. It remains to explain the \kappa. B cells produce antibodies each of which consists of two symmetrical halves. Each half consists of a light and a heavy chain. The light chain comes in two variants called \kappa and \lambda. The choice which of these a cell uses seems to be fairly random. The work in the Baltimore lab had found out that NF\kappaB could skew the ratio. I found a video by Baltimore from 2001 about NF\kappaB. This is probably quite out of date by now but it contained one thing which I found interesting. Under certain circumstances it can happen that a constant stimulus causing activation of NF\kappaB leads to oscillations in the concentration. In the video the speaker mentions ‘odd oscillations’ and comments ‘but that’s for mathematicians to enjoy themselves’. It seems that he did not believe these oscillations to be biologically important. There are reasons to believe that they might be important and I will try to explain why. At the very least it will allow me to enjoy myself.

Let me explain the usual story about how NF\kappaB is activated. There are lots of animated videos on Youtube illustrating this but I prefer a description in words. Normally NF\kappaB is found in the cytosol bound to an inhibitor I\kappaB. Under certain circumstances a complex of proteins called IKK forms. The last K stands for kinase and IKK phosphorylates I\kappaB. This causes I\kappaB to be ubiquinated and thus marked for degradation (cf. the discussion of ubiquitin here). When it has been destroyed NF\kappaB is liberated, moves to the nucleus and binds to DNA. What are the circumstances mentioned above? There are many alternatives. For instance TNF\alpha binds to its receptor, or something stimulates a toll-like receptor. The details are not important here. What is important is that there are many different signals which can lead to the activation of NF\kappaB. What genes does NF\kappaB bind to when it is activated? Here again there are many possibilities. Thus there is a kind of bow tie configuration where there are many inputs and many outputs which are connected to a single channel of communication. So how is it possible to arrange that when one input is applied, e.g. TNF\alpha the right genes are switched on while another input activates other genes through the same mediator NF\kappaB? One possibility is cross-talk, i.e. that this signalling pathway interacts with others. If this cannot account for all the specificity then the remaining possibility is that information is encoded in the signal passing through NF\kappaB itself. For example, one stimulus could produce a constant response while another causes an oscillatory one. Or two stimuli could cause oscillatory responses with different frequencies. Evidently the presence of oscillations in the concentration of NF\kappaB presents an opportunity for encoding more information than would otherwise be possible. To what extent this really happens is something where I do not have an overview at the moment. I want to learn more. In any case, oscillations have been observed in the NF\kappaB system. The primary thing which has been observed to oscillate is the concentration of NF\kappaB in the nucleus. This oscillation is a consequence of the movement of the protein between the cytosol and the nucleus. There are various mathematical models for describing these oscillations. As usual in modelling phenomena in cell biology there are models which are very big and complicated. I find it particularly interesting when some of the observations can be explained by a simple model. This is the case for NF\kappaB where a three-dimensional model and an explanation of its relations to the more complicated models can be found in a paper by Krishna, Jensen and Sneppen (PNAS 103, 10840). In the three-dimensional model the unknowns are the concentrations of NF\kappaB in the nucleus, I\kappaB in the cytoplasm and mRNA coding for I\kappaB. The oscillations in normal cells are damped but sustained oscillations can be seen in mutated cells or corresponding models.

What is the function of NF\kappaB? The short answer is that it has many. On a broad level of description it plays a central role in the phenomenon of inflammation. In particular it leads to production of the cytokine IL-17 which in turn, among other things, stimulates the production of anti-microbial peptides. When these things are absent it leads to a serious immunodeficiency. In one variant of this there is a mutation in the gene coding for NEMO, which is one of the proteins making up IKK. A complete absence of NEMO is fatal before birth but people with a less severe mutation in the gene do occur. There are symptoms due to things which took place during the development of the embryo and also immunological problems, such as the inability to deal with certain bacteria. The gene for NEMO is on the X chromosome so that this disease is usually limited to boys. More details can be found in the book of Geha and Notarangelo mentioned in  a previous post.

Stability of steady states in models of the Calvin cycle

April 25, 2016

I have just written a paper with Stefan Disselnkötter on stationary solutions of models for the Calvin cycle and their stability. There we concentrate on the simplest models for this biological system. There were already some analytical results available on the number of positive stationary solutions (let us call them steady states for short), with the result that this number is zero, one or two in various circumstances. We were able to extend these results, in particular showing that in a model of Zhu et. al. there can be two steady states or, in exceptional cases, a continuum of steady states. This is at first sight surprising since those authors stated that there is at most one steady state. However they impose the condition that the steady states should be ‘physiologically feasible’. In fact for their investigations, which are done by means of computer calculations, they assume among other things that certain Michaelis constants which occur as parameters in the system have specific numerical values. This assumption is biologically motivated but at the moment I do not understand how the numbers they give follow from the references they quote. In any case, if these values are assumed our work gives an analytical proof that there is at most one steady state.

While there are quite a lot of results in the literature on the number of steady states in systems of ODE modelling biochemical systems there is much less on the question of the stability of these steady states. It was a central motivation of our work to make some progress in this direction for the specific models of the Calvin cycle and to develop some ideas to approaching this type of question more generally. One key idea is that if it can be shown that there is bifurcation with a one-dimensional centre manifold this can be very helpful in getting information on the stability of steady states which arise in the bifurcation. Given enough information on a sufficient number of derivatives at the bifurcation point this is a standard fact. What is interesting and perhaps less well known is that it may be possible to get conclusions without having such detailed control. One type of situation occurring in our paper is one where a stable solution and a saddle arise. This is roughly the situation of a fold bifurcation but we do not prove that it is generic. Doing so would presumably involve heavy calculations.

The centre manifold calculation only controls one eigenvalue and the other important input in order to see that there is a stable steady state for at least some choice of the parameters is to prove that the remaining eigenvalues have negative real parts. This is done by considering a limiting case where the linearization simplifies and then choosing parameters close to those of the limiting case. The arguments in this paper show how wise it can be to work with the rates of the reactions as long as possible, without using species concentrations. This kind of approach is popular with many people – it has just taken me a long time to get the point.

The advanced deficiency algorithm

January 23, 2016

Here I discuss another tool for analysing chemical reaction networks of deficiency greater than one. This is the Advanced Deficiency Algorithm developed by Feinberg and Ellison. It seems that the only direct reference for the mathematical proofs is Ellison’s PhD thesis. There is a later PhD thesis by Haixia Ji in which she introduces an extension of this called the Higher Deficiency Algorithm and where some of the ideas of Ellison are also recapitulated. In my lecture course, which ends next week, I will only have time to discuss the structure of the algorithm and give an extended example without proving much.

The Advanced Deficiency Algorithm has a general structure which is similar to that of the Deficiency One Algorithm. In some cases it can rule out multistationarity. Otherwise it gives rise to several sets of inequalities. If one of these has a solution then there is multistationarity and if none of them does there is no multistationarity. It is not clear to me if this is really an algorithm which is guaranteed to give a diagostic test in all cases. I think that this is probably not the case and that one of the themes of Ji’s thesis is trying to improve on this. An important feature of this algorithm is that the inequalities it produces are in general nonlinear and thus may be much more difficult to analyse than the linear inequalities obtained in the case of the Deficiency One Algorithm.

Now I have come to the end of my survey of deficiency theory for chemical reaction networks. I feel I have learned a lot and now is the time to profit from that by applying these techniques. The obvious next step is to try out the techniques on some of my favourite biological examples. Even if the result is only that I see why the techniques do not give anything interesting in this cases it will be useful to understand why. Of course I hope that I will also find some positive results.

Elementary flux modes

January 16, 2016

Elementary flux modes are a tool which can be used to take a reaction network of deficiency greater than one and produce subnetworks of deficiency one. If it can be shown that one of these subnetworks admits multistationarity then it can sometimes be concluded that the original network also does so. There are two conditions that need to be checked in order to allow this conclusion. The first is that the network under consideration must satisfy the condition t=l. The second is that genericity conditions must be satisfied which consist of requiring the non-vanishing of the determinants of certain matrices. If none of the subnetworks admit multistationarity then no conclusion is obtained for the full network. The core element of the proof is an application of the implicit function theorem. These conclusions are contained in a paper of Conradi et. al. (Proc. Nat. Acad. Sci. USA 104, 19175).

One of the ways of writing the condition for stationary solutions is Nv(c)=0, where as usual N is the stoichiometric matrix and v(c) is the vector of reaction rates. Since v(c) is positive this means that we are looking for a positive element of the kernel of N. This suggests that it is interesting to look at the cone K which is the intersection of the kernel of N with the non-negative orthant. According to a general theorem of convex analysis K consists of the linear combinations with non-negative coefficients of a finite set of vectors which have a mininum (non-zero) number of non-zero components. In the case of reaction networks these are the elementary flux modes. Recalling that N=YI_a we see that positive vectors in the kernel of the incidence matrix I_a are a special type of elementary flux modes. Those which are not in the kernel of I_a are called stoichiometric generators. Each stoichiometric generator defines a subnetwork where those reaction constants of the full network are set to zero where the corresponding reactions are not in the support of the generator. It is these subnetworks which are the ones mentioned above in the context of multistationarity. The application of the implicit function theorem involves using a linear transformation to introduce adapted coordinates. Roughly speaking the new coordinates are of three types. The first are conserved quantities for the full network. The second are additional conserved quantities for the subnetwork, complementing those of the full network. Finally the third type represents quantities which are dynamical even for the subnetwork.

Here are some simple examples. In the extended Michaelis-Menten description of a single reaction there is just one elementary flux mode (up to multiplication by a positive constant) and it is not a stoichiometric generator. In the case of the simple futile cycle there are three elementary flux modes. Two of these, which are not stoichoimetric generators correspond to the binding and unbinding of one of the enzymes with its substrate. The third is a stoichoimetric generator and the associated subnetwork is obtained by removing the reactions where a substrate-enzyme complex dissociates into its original components. The dual futile cycle has four elementary flux modes of which two are stoichiometric generators. In the latter case we get the (stoichiometric generators of the) two simple futile cycles contained in the network. Of course these are not helpful for proving multistationarity. Another type of example is given by the simple models for the Calvin cycle which I have discussed elsewhere. The MA system considered there has two stoichiometric generators with components (3,6,6,1,3,0,1) and (3,5,5,1,3,1,0). I got these by working backwards from corresponding modes for the MM-MA system found by Grimbs et. al. This is purely on the level of a vague analogy. I wish I had a better understanding of how to get this type of answer more directly. Those authors used those modes to prove bistability for the MM-MA system so that this is an example where this machinery produces real results.

The deficiency one algorithm

January 7, 2016

Here I continue the discussion of Chemical Reaction Network Theory with the Deficiency One Algorithm. As its name suggests this is an algorithm for investigating whether a reaction network allows multistationarity or not. It is also associated to a theorem. This theorem says that the satisfaction of certain inequalities associated to a reaction network is equivalent to the existence of a choice of reaction constants for which the corresponding system of equations with mass action kinetics has two distinct positive stationary solutions in a stoichiometric compatibility class. This is true in the context of networks of deficiency one which satisfy some extra conditions. The property of satisfying these, called regularity, seems to be relatively weak. In fact, as in the case of the Deficiency One Theorem, there is a second related result where a pair of zeroes of the right hand side f in a stoichiometric compatibility class is replaced by a vector in the kernel of Df which is tangent to that class.

An example where this theory can be applied is double phosphorylation in the processive case. The paper of Conradi et. al. cited in the last post contains the statement that in this system the Deficiency One Algorithm implies that multistationarity is not possible. For this it refers to the Chemical Reaction Network Toolbox, a software package which implements the calculations of the theorem. In my course I decided for illustrative purposes to carry out these calculations by hand. It turned out not to be very hard. The conclusion is that multistationarity is not possible for this system but the general machinery does not address the question of whether there are any positive stationary solutions. I showed that this is the case by checking by hand that \omega-limit points on the boundary of the positive orthant are impossible. The conclusion then follows by a well-known argument based on the Brouwer fixed point theorem. This little bit of practical experience with the Deficiency One Algorithm gave me the impression that it is really a powerful tool. At this point it is interesting to note a cross connection to another subject which I have discussed in this blog. It is a model for the Calvin cycle introduced by Grimbs et. al. These authors noted that the Deficiency One Algorithm can be applied to this system to show that it does not allow multistationarity. They do not present the explicit calculations but I found that they are not difficult to do. In this case the equations for stationary solutions can be solved explicitly so that using this tool is a bit of an overkill. It neverless shows the technique at work in another example.

Regularity consists of three conditions. The first (R1) is a necessary condition that there be any positive solutions at all. If it fails we get a strong conclusion. The second (R2) is the condition l=t familiar from the Deficiency One Theorem. (R3) is a purely graph theoretic condition on the network. A weakly reversible network which satisfies (R3) is reversible. Reading this backwards, if a network is weakly reversible but not reversible then the theorem does not apply. The inequalities in the theorem depend on certain partitions of a certain set of complexes (the reactive complexes) into three subsets. What is important is whether the inequalities hold for all partitions of a network or whether there is at least one partition for which they do not hold. The proof of the theorem uses a special basis of the kernel of YA_\kappa where \kappa is a function on \cal C constructed from the reaction constants. In the context of the theorem this space has dimension t+1 and t of the basis vectors come from a special basis of A_\kappa of the type which already comes up in the proof of the Deficiency Zero Theorem.

An obvious restriction on the applicability of the Deficiency One Algorithm is that it is only applicable to networks of deficiency one. What can be done with networks of higher deficiency? One alternative is the Advanced Deficiency Algorithm, which is implemented in the Chemical Reaction Network Toolbox. A complaint about this method which I have seen several times in the literature is that it is not able to treat large systems – apparently the algorithm becomes unmanageable. Another alternative uses the notion of elementary flux modes which is the next topic I will cover in my course. It is a way of producing certain subnetworks of deficiency one from a given network of higher deficiency. The subnetworks satisfy all the conditions needed to apply the Deficiency One Algorithm except perhaps t=l.

The deficiency one theorem

December 2, 2015

Here I continue with the discussion of chemical reaction network theory begun in the previous post. After having presented a proof of the Deficiency Zero Theorem in my course I proceeded to the Deficiency One Theorem following the paper of Feinberg in Arch. Rat. Mech. Anal. 132, 311. If we have a reaction network we can consider each linkage class as a network in its own right and thus we can define its deficiency. If we denote the deficiency of the full network by \delta and the deficiency of the ith linkage class by \delta_i then in general \delta\ge\sum_i \delta_i. The first hypothesis of the Deficiency One Theorem is that the deficiencies of the linkage classes are no greater than one. The second is that equality holds in the inequality relating the deficiency of the network to those of its linkage classes. The stoichiometric subspace S of the full network is the sum of the corresponding spaces S_i for the linkage classes. The sum is direct precisely when equality holds in the inequality for the deficiencies. The third condition is that there is precisely one terminal strong linkage class in each linkage class (t=l). The first conclusion is that if there exists a positive stationary solution there is precisely one stationary solution in each stoichiometric compatibility class. The second is that if the network is weakly reversible there exists a positive stationary solution. Thus two of the conclusions of the Deficiency Zero Theorem have direct analogues in this case. Others do not. There is no statement about the stability of the stationary solutions. Networks which are not weakly reversible but satisfy the three conditions of the Deficiency One Theorem may have positive stationary solutions. In the paper of Feinberg the proof of the Deficiency One Theorem is intertwined with that of another result. It says that the linearization about a stationary solution of the restriction of the system to a stoichiometric compatibility class has trivial kernel. A related result proved in the same paper says that for a weakly reversible network of deficiency zero each stationary solution is a hyperbolic sink within its stoichiometric class. Statements of this type ensure that these stationary solutions possess a certain structural stability.

In the proof of the Deficiency One Theorem the first step (Step 1) is to show that when the assumptions of the theorem are satisfied and there is positive stationary solution c^* then the set of stationary solutions is equal to the set of points for which \log c-\log c^* lies in the orthogonal complement of the stoichiometric subspace. From this the conclusion that there is exactly one stationary solution in each stoichiometric compatibility class follows just as in the proof of the Deficiency Zero Theorem (Step 2). To complete the proof it then suffices to prove the existence of a positive stationary solution in the weakly reversible case (Step 3). In Step 1 the proof is reduced to the case where the network has only one linkage class by regarding the linkage classes of the original network as networks in their own right. In this context the concept of a partition of a network (\cal S,\cal C,\cal R) is introduced. This is a set of subnetworks ({\cal S},{\cal C}^i,{\cal R}^i). The set of species is unchanged. The set of reactions \cal R is a disjoint union of the {\cal R}^i) and {\cal C}^i is the set of complexes occurring in the reactions contained in {\cal R}^i). The partition is called direct if the stoichiometric subspace of the full network is the direct sum of those of the subnetworks. The most important example of a partition of a network in the present context is that given by the linkage classes of any network. That it is direct is the second condition of the Deficiency One Theorem. The other part of Step 1 of the proof is to show that the statement holds in the case that there is only one linkage class. The deficiency is then either zero or one. Since the case of deficiency zero is already taken care of by the Deficiency Zero Theorem we can concentrate on the case where \delta=1. Then the dimension of {\rm ker} A_k is one and that of {\rm ker} (YA_k) is two. The rest of Step 1 consists of a somewhat intricate algebraic calculation in this two-dimensional space. It remains to discuss Step 3. In this step the partition given by the linkage classes is again used to reduce the problem to the case where there is only one linkage class. The weak reversibility is preserved by this reduction. Again we can assume without loss of generality that \delta=1. The subspace U={\rm im} Y^T+{\rm span}\omega_{\cal C} is a hyperplane in F({\cal C}). We define \Gamma to be the set of functions of the form \log a with a a positive element of {\rm ker} (YA_k). The desired stationary solution is obtained as a point of the intersection of U with \Gamma. To show that this intersection is non-empty it is proved that there are points of \Gamma on both sides of U. This is done by a careful examination of the cone of positive elements of {\rm ker} (YA_k).

To my knowledge the number of uses of the Deficiency Zero Theorem and the Deficiency One Theorem in problems coming from applications is small. If anyone reading this has other information I would like to hear it. I will now list the relevant examples I know. The Deficiency Zero Theorem was applied by Sontag to prove asymptotic stability for McKeithan’s kinetic proofreading model of T cell activation. I applied it to the model of Salazar and Höfer for the dynamics of the transcription factor NFAT. Another potential application would be the multiple futile cycle. This network is not weakly reversible. The simple futile cycle has deficiency one and the dual futile cycle deficiency two. The linkage classes have deficiency zero in both cases. Thus while the first and third conditions of the Deficiency One Theorem are satisfied the second is not. Replacing the distributive phosphorylation in the multiple futile cycle by processive phosphorylation may simplify things. This has been discussed by Conradi et. al., IEE Proc. Systems Biol. 152, 243. In the case of two phosphorylation steps the system obtained is of deficiency one but the linkage classes have deficiency zero so that condition (ii) is still not satisfied. It seems that the Deficiency One Algorithm may be more helpful and that will be the next subject I cover in my course.

Feinberg’s proof of the deficiency zero theorem

November 23, 2015

I discussed the deficiency zero theorem of chemical reaction network theory (CRNT) in a previous post. (Some further comments on this can be found here and here.) This semester I am giving a lecture course on chemical reaction network theory. Lecture notes are growing with the course and can be found in German and English versions on the web page of the course. The English version can also be found here. Apart from introductory material the first main part of the course was a proof of the Deficiency Zero Theorem. There are different related proofs in the literature and I have followed the approach in the classic lecture notes of Feinberg on the subject closely. The proof in those notes is essentially self-contained apart from one major input from a paper of Feinberg and Horn (Arch. Rat. Mech. Anal. 66, 83). In this post I want to give a high-level overview of the proof.

The starting point of CRNT is a reaction network. It can be represented by a directed graph where the nodes are the complexes (left or right hand sides of reactions) and the directed edges correspond to the reactions themselves. The connected components of this graph are called the linkage classes of the network and their number is usually denoted by l. If two nodes can be connected by oriented paths in both directions they are said to be strongly equivalent. The corresponding equivalence classes are called strong linkage classes. A strong linkage class is called terminal if there is no directed edge leaving it. The number of terminal strong linkage classes is usually denoted by t. From the starting point of the network making the assumption of mass action kinetics allows a system of ODE \dot c=f(c) to be obtained in an algorithmic way. The quantity c(t) is a vector of concentrations as a function of time. Basic mathematical objects involved in the definition of the network are the set \cal S of chemical species, the set \cal C of complexes and the set \cal R of reactions. An important role is also played by the vector spaces of real-valued functions on these finite sets which I will denote by F({\cal S}), F({\cal C}) and F({\cal R}), respectively. Using natural bases they can be identified with R^m, R^n and R^r. The vector c(t) is an element of F({\cal S}). The mapping f from F({\cal S}) to itself can be written as a composition of three mappings, two of them linear, f=YA_k\Psi. Here Y, the complex matrix, is a linear mapping from F({\cal C}) to F({\cal S}). A_k is a linear mapping from F({\cal C}) to itself. The subscript k is there because this matrix is dependent on the reaction constants, which are typically denoted by k. It is also possible to write f in the form Nv where v describes the reaction rates and N is the stoichiometric matrix. The image of N is called the stoichiometric subspace and its dimension, the rank of the network, is usually denoted by s. The additive cosets of the stoichiometric subspace are called stoichiometric compatibility classes and are clearly invariant under the time evolution. Finally, \Psi is a nonlinear mapping from F({\cal S}) to F({\cal C}). The mapping \Psi is a generalized polynomial mapping in the sense that its components are products of powers of the components of c. This means that \log\Psi depends linearly on the logarithms of the components of c. The condition for a stationary solution can be written as \Psi(c)\in {\rm ker} (YA_k). The image of \Psi is got by exponentiating the image of a linear mapping. The matrix of this linear mapping in natural bases is Y^T. Thus in looking for stationary solutions we are interested in finding the intersection of the manifold which is the image of \Psi with the kernel of YA_k. The simplest way to define the deficiency of the network is to declare it to be \delta=n-l-s. A fact which is not evident from this definition is that \delta is always non-negative. In fact \delta is the dimension of the vector space {\rm ker} Y\cap {\rm span}(\Delta) where \Delta is the set of complexes of the network. An alternative concept of deficiency, which can be found in lecture notes of Gunawardena, is the dimension \delta' of the space {\rm ker} Y\cap{\rm im} A_k. Since this vector space is a subspace of the other we have the inequality \delta'\le\delta. The two spaces are equal precisely when each linkage class contains exactly one terminal strong linkage class. This is, in particular, true for weakly reversible networks. The distinction between the two definitions is often not mentioned since they are equal for most networks usually considered.

If c^* is a stationary solution then A_k\Psi (c^*) belongs to {\rm ker} Y\cap{\rm im} A_k. If \delta'=0 (and in particular if \delta=0) then this means that A_k\Psi (c^*)=0. In other words \Psi (c^*) belongs to the kernel of A_k. Stationary solutions of this type are called complex balanced. It turns out that if c^* is a complex balanced stationary solution the stationary solutions are precisely those points c for which \log c-\log c_* lies in the orthogonal complement of the stoichiometric subspace. It follows that whenever we have one solution we get a whole manifold of them of dimension n-s. It can be shown that each manifold of this type meets each stoichiometric class in precisely one point. This is proved using a variational argument and a little convex analysis.

It is clear from what has been said up to now that it is important to understand the positive elements of the kernel of A_k. This kernel has dimension t and a basis each of whose elements is positive on a terminal strong linkage class and zero otherwise. Weak reversibility is equivalent to the condition that the union of the terminal strong linkage classes is the set of all complexes. It can be concluded that when the network is not weakly reversible there exists no positive element of the kernel of A_k. Thus for a network which is not weakly reversible and has deficiency zero there exist no positive stationary solutions. This is part of the Deficiency Zero Theorem. Now consider the weakly reversible case. There a key statement of the Deficiency Zero Theorem is that there exists a complex balanced stationary solution c^*. Where does this c^* come from? We sum the vectors in the basis of {\rm ker} A_k and due to weak reversibility this gives something which is positive. Then we take the logarithm of the result. When \delta=0 this can be represented as a sum of two contributions where one is of the form Y^T z. Then c^*=e^z. A further part of the deficiency zero theorem is that the stationary solution c^* in the weakly reversible case is asymptotically stable. This is proved using the fact that for a complex balanced stationary solution the function h(c)=\sum_{s\in\cal S}[c(s)(\log c(s)-\log c^*(s)-1)+c(s^*)] is a Lyapunov function which vanishes for c=c^*

Siphons in reaction networks

October 8, 2015

The concept of a siphon is one which I have been more or less aware of for quite a long time. Unfortunately I never had the impression that I had understood it completely. Given the fact that it came up a lot in discussions I was involved in and talks I heard last week I thought that the time had come to make the effort to do so. It is of relevance for demonstrating the property of persistence in reaction networks. This is the property that the \omega-limit points of a positive solution are themselves positive. For a bounded solution this is the same as saying that the infima of all concentrations at late times are positive. The most helpful reference I have found for these topics is a paper of Angeli, de Leenheer and Sontag in a proceedings volume edited by Queinnec et. al.

There are two ways of formulating the definition of a siphon. The first is more algebraic, the second more geometric. In the first the siphon is defined to be a set Z of species with the property that whenever one of the species in Z occurs on the right hand side of a reaction one of the species in Z occurs on the left hand side. Geometrically we replace Z by the set L_Z of points of the non-negative orthant which are common zeroes of the elements of Z, thought of as linear functions on the species space. The defining property of a siphon is that L_Z is invariant under the (forward in time) flow of the dynamical system describing the evolution of the concentrations. Another way of looking at the situation is as follows. Consider a point of L_Z. The right hand side of the evolution equations of one of the concentrations belonging to Z is a sum of positive and negative terms. The negative terms automatically vanish on L_Z and the siphon condition is what is needed to ensure that the positive terms also vanish there. Sometimes minimal siphons are considered. It is important to realize that in this case Z is minimal. Correspondingly L_Z is maximal. The convention is that the empty set is excluded as a choice for Z and correspondingly the whole non-negative orthant as a choice for L_Z. What is allowed is to choose Z to be the whole of the species space which means that L_Z is the origin. Of course whether this choice actually defines a siphon depends on the particular dynamical system being considered.

If x_* is an \omega-limit point of a positive solution but is not itself positive then the set of concentrations which are zero at that point is a siphon. In particular stationary solutions on the boundary are contained in siphons. It is remarked by Shiu and Sturmfels (Bull. Math. Biol. 72, 1448) that for a network with only one linkage class if a siphon contains one stationary solution it consists entirely of stationary solutions. To see this let x_* be a stationary solution in the siphon Z. There must be some complex y belonging to the network which contains an element of Z. If y' is another complex then there is a directed path from y' to y. We can follow this path backwards from y and conclude successively that each complex encountered contains an element of Z. Thus y' contains an element of Z and since y' was arbitrary all complexes have this property. This means that all complexes vanish at x_* so that x_* is a stationary solution.

Siphons can sometimes be used to prove persistence. Suppose that Z is a siphon for a certain network so that the points of Z are potential \omega-limit points of solutions of the ODE system corresponding to this network. Suppose further that A is a conserved quantity for the system which is a linear combination of the coordinates with positive coefficents. For a positive solution the quantity A has a positive constant value along the solution and hence also has the same value at any of its \omega-limit points. It follows that if A vanishes on Z then no \omega-limit point of that solution belongs to Z. If it is possible to find a conserved quantity A of this type for each siphon of a given system (possibly different conserved quantities for different siphons) then persistence is proved. For example this strategy is used in the paper of Angeli et al. to prove persistence for the dual futile cycle. The concept of persistence is an important one when thinking about the general properties of reaction networks. The persistence conjecture says that any weakly reversible reaction network with mass action kinetics is persistent (possibly with the additional assumption that all solutions are bounded). In his talk last week Craciun mentioned that he is working on proving this conjecture. If true it implies the global attractor conjecture. It also implies a statement claimed in a preprint of Deng et. al. (arXiv:1111.2386) that a weakly reversible network has a positive stationary solution in any stoichiometric compatobility class. This result has never been published and there seems to be some doubt as to whether the proof is correct.

Trip to the US

October 5, 2015

Last week I visited a few places in the US. My first stop was Morgantown, West Virginia where my host was Casian Pantea. There I had a lot of discussions with Casian and Carsten Conradi on chemical reaction network theory. This synergized well with the work I have recently been doing preparing a lecture course on that subject which I will be giving in the next semester. I gave a talk on MAPK and got some feedback on that. It rained a lot and there was not much opportunity to do anything except work. One day on the way to dinner while it was relatively dry I saw a Cardinal and I fortunately did have my binoculars with me. On Wednesday afternoon I travelled to New Brunswick and spent most of Thursday talking to Eduardo Sontag at Rutgers. It was a great pleasure to talk to an excellent mathematician who also knows a lot about immunology. He and I have a lot of common interests which is in part due to the fact that I was inspired by several of his papers during the time I was getting into mathematical biology. I also had the opportunity to meet Evgeni Nikolaev who told me a variety of interesting things. They concerned bifurcation theory in general, its applications to the kinds of biological models I am interested in and his successes in applying mathematical models to understanding concrete problems in biomedical research such as the processes taking place in tuberculosis. My personal dream is to see a real coming together of mathematics and immunology and that I have the chance to make a contribution to that process.

On Friday I flew to Chicago in order to attend an AMS sectional meeting. I had been in Chicago once before but that is many years ago now. I do remember being impressed by how much Lake Michigan looks like the sea, I suppose due to the structure of the waves. This impression was even stronger this time since there were strong winds whipping up the waves. Loyola University, the site of the meeting, is right beside the lake and it felt like home for me due to the combination of wind, waves and gulls. The majority of those were Ring-Billed Gulls which made it clear which side of the Atlantic I was on. There were also some Herring Gulls and although they might have been split from those on the other side of the Atlantic by the taxonomists I did not notice any difference. It was the first time I had been at an AMS sectional meeting and my impression was that the parallel sessions were very parallel, in other words in no danger of meeting. Most of the people in our session were people I knew from the conferences I attended in Charlotte and in Copenhagen although I did make a couple of new acquaintances, improving my coverage of the reaction network community.

In a previous post I mentioned Gheorghe Craciun’s ideas about giving the deficiency of a reaction network a geometric interpretation, following a talk of his in Copenhagen. Although I asked him questions about this on that occasion I did not completely understand the idea. Correspondingly my discussion of the point here in my blog was quite incomplete. Now I talked to him again and I believe I have finally got the point. Consider first a network with a single linkage class. The complexes of the network define points in the species space whose coordinates are the stoichiometric coefficients. The reactions define oriented segments joining the educt complex to the product complex of each reaction. The stoichiometric subspace is the vector space spanned by the differences of the complexes. It can also be considered as a translate of the affine subspace spanned by the complexes themselves. This makes it clear that its dimension s is at most n-1, where n is the number of complexes. The number s is the rank of the stoichiometric matrix. The deficiency is n-1-s. At the same time s\le m. If there are several linkage classes then the whole space has dimension at most n-l, where l is the number of linkage classes. The deficiency is n-l-s. If the spaces corresponding to the individual linkage classes have the maximal dimension allowed by the number of complexes in that class and these spaces are linearly independent then the deficiency is zero. Thus we see that the deficiency is the extent to which the complexes fail to be in general position. If the species and the number of complexes have been fixed then deficiency zero is seen to be a generic condition. On the other hand fixing the species and adding more complexes will destroy the deficiency zero condition since then we are in the case n-l>m so that the possibility of general position is excluded. The advantage of having this geometric picture is that it can often be used to read off the deficiency directly from the network. It might also be used to aid in constructing networks with a desired deficiency.