## Archive for the ‘dynamical systems’ Category

### Conference on reaction networks and population dynamics in Oberwolfach

July 9, 2017

This is a belated report on a conference in Oberwolfach I attended a couple of weeks ago. The title includes two elements. The first held no suprises for me but the second was rather different from what I had expected. My expectation was that it would be about the evolution of populations of organisms. In fact it was rather focussed on models related to genetics, in other words with the question of how certain genetic traits spread through a population.

There was one talk which did have a connection to population biology in way closer to what I had expected. It happens all the time that ecosystems are damaged by exotic species imported, deliberately or by accident, from other parts of the world. There are also well-known stories of the type that to try to control exotic species number one exotic species number two is introduced and is itself very harmful. It is nice to hear an example where this kind of introduction of an exotic species was very successful. It is the case of the cassava plant which was introduced from South America to Africa and became a staple food there. Then an insect from South America (species number one) called the mealy bug was introduced accidentally and caused enormous damage. Finally an ecologist called Hans Herren introduced a parasitic wasp (species number two) from South America, restoring the food supply and saving numerous lives (often the number 20 million is quoted). More details of this story can be found here.

I want to mention one statement made in the talk of Gheorghe Craciun in Oberwolfach which I found intriguing. I might have heard this before but it did not stick in my mind properly. The statement is that the set of dynamical systems which possess a complex balanced steady state is a variety of codimension $\delta$, where $\delta$ is the deficiency. There seemed to be some belief in the audience that this variety is actually a smooth manifold. On one afternoon we had something similar to the breakout sessions in Banff. I suggested the topic for one of these, which was Lyapunov functions. The idea was to compare classes of Lyapunov functions which people working on different classes of dynamical systems knew. This certainly did not lead to any breakthrough but I think it did lead to a useful exchange of information. I documented the discussion for my own use and I think I could profit by following some of the leads there.

To finish I want to mention a claim made by Ankit Gupta in his talk. It did not sound very plausible to me but I expect that it at least contains a grain of truth. He said that these days more papers are published on $NF\kappa B$ than on all of mathematics.

### Conference on mathematical analysis of biological interaction networks at BIRS

June 9, 2017

I have previously written a post concerning a meeting at the Banff International Research Station (BIRS). This week I am at BIRS again. Among the topics of the talks were stochastic chemical reaction networks, using reaction networks in cells as computers and the area of most direct relevance to me, multiple steady states and their stability in deterministic CRN. Among the most popular examples occurring in the talks in the latter area were the multiple futile cycle, the MAPK cascade and the EnvZ/OmpR system. In addition to the talks there was a type of event which I had never experienced before called breakout sessions. There the participants split into groups to discuss different topics. The group I joined was concerned with oscillations in phosphorylation cycles.

In the standard dual futile cycle we have a substrate which can be phosphorylated up to two times by a kinase and dephosphorylated again by a phosphatase. It is assumed that the (de-)phosphorylation is distributive (the number of phosphate groups changes by one each time a substrate binds to an enzyme) and sequential (the phosphate groups are added in one order and removed in the reverse order). A well-known alternative to this is processive (de-)phosphorylation where the number of phosphate groups changes by two in one encounter between a substrate and an enzyme. It is known that the double phosphorylation system with distributive and sequential phosphorylation admits reaction constants for which there are three steady states, two of which are stable. (From now on I only consider sequential phosphorylation here.) By contrast the corresponding system with processive phophorylation always has a unique steady state. Through the talk of Anne Shiu here I became aware of the following facts. In a paper by Suwanmajo and Krishnan (J. R. Soc. Interface 12:20141405) it is stated that in a mixed model with distributive phosphorylation and processive dephosphorylation periodic solutions occur as a result of a Hopf bifurcation. The paper does not present an analytical proof of this assertion.

It is a well-known open question, whether there are periodic solutions in the case that the modificiations are all distributive. It has been claimed in a paper of Errami et. al. (J. Comp. Phys. 291, 279) that a Hopf bifurcation had been discovered in this system but the claim seems to be unjustified. In our breakout sessions we looked at whether oscillations might be exported from the mixed model to the purely distributive model. We did not get any definitive results yet. There were also discussions on effective ways of detecting Hopf bifurcations, for instance by using Hurwitz determinants. It is well-known that oscillations in the purely distributive model, if they exist, do not persist in the Michaelis-Menten limit. I learned from Anne Shiu that it is similarly the case that the oscillations in the mixed model are absent from the Michaelis-Menten system. This result came out of some undergraduate research she supervised. Apart from these specific things I learned a lot just from being in the environment of these CRN people.

Yesterday was a free afternoon and I went out to look for some birds. I saw a few things which were of interest to me, one of which was a singing Tennessee warbler. This species has a special significance for me for the following reason. Many years ago when I still lived in Orkney I got an early-morning phone call from Eric Meek, the RSPB representative. He regularly checked a walled garden at Graemeshall for migrants. On that day he believed he had found a rarity and wanted my help in identifying it and, if possible, catching it. We did catch it and it turned out to be a Tennessee warbler, the third ever recorded in Britain. That was big excitement for us. I had not seen Eric for many years and I was sad to learn now that he died a few months ago at a relatively young age. The name of this bird misled me into thinking that it was at home in the southern US. In fact the name just came from the fact that the first one to be described was found in Tennessee, a migrant. The breeding range is much further north, especially in Canada. Thus it is quite appropriate that I should meet it here.

### The Routh-Hurwitz criterion

May 7, 2017

I have been aware of the Routh-Hurwitz criterion for stability for a long time and I have applied it in three dimensions in my research and tried to apply it in four. Unfortunately I never felt that I really understood it completely. Here I want to finally clear this up. A source which I found more helpful than other things I have seen is https://www.math24.net/routh-hurwitz-criterion/. One problem I have had is that the Hurwitz matrices, which play a central role in this business, are often written in a form with lots of … and I was never sure that I completely understood the definition. I prefer to have a definite algorithm for constructing these matrices. The background is that we would like to understand the stability of steady states of a system of ODE. Suppose we have a system $\dot x=f(x)$ and a steady state $x_0$, i.e. a solution of $f(x_0)=0$. It is well-known that this steady state is asymptotically stable if all eigenvalues $\lambda$ of the linearization $A=Df(x_0)$ have negative real parts. This property of the eigenvalues is of course a property of the roots of the characteristic equation $\det(A-\lambda I)=a_0\lambda^n+\ldots+a_{n-1}\lambda+a_n=0$. It is always the case here that $a_0=1$ but I prefer to deal with a general polynomial with real coefficients $a_i, 0\le i\le n$ and a criterion for the situation where all its roots have negative real parts. It is tempting to number the coefficients in the opposite direction, so that, for instance, $a_n$ becomes $a_0$ but I will stick to this convention. Note that it is permissible to replace $a_k$ by $a_{n-k}$ in any criterion of this type since if we multiply the polynomial by $\lambda^{-n}$ we get a polynomial in $\lambda^{-1}$ where the order of the coefficients has been reversed. Moreover, if the real part of $\lambda$ is non-zero then it has the same sign as the real part of $\lambda^{-1}$. I find it important to point this out since different authors use different conventions for this. It is convenient to formally extend the definition of the $a_i$ to the integers so that these coefficients are zero for $i<0$ and $i>n$.

For a fixed value of $n$ the Hurwitz matrix is an $n$ by $n$ matrix defined as follows. The $j$th diagonal element is $a_j$, with $1\le j\le n$. Starting from a diagonal element and proceeding to the left along a row the index increases by one in each step. Similarly, proceeding to the right along a row the index decreases by one. In the ranges where the index is negative or greater than $n$ the element $a_n$ can be replaced by zero. The leading principal minors of the Hurwitz matrix, in other words the determinants of the submatrices which are the upper left hand corner of the original matrix, are the Hurwitz determinants $\Delta_k$. The Hurwitz criterion says that the real parts of all roots of the polynomial are negative if and only if $a_0>0$ and $\Delta_k>0$ for all $1\le k\le n$. Note that a necessary condition for all roots to have negative real parts is that all $a_i$ are positive. Now $\Delta_n=a_n\Delta_{n-1}$ and so the last condition can be replaced by $a_n>0$. Note that the form of the $\Delta_k$ does not depend on $n$. For $n=2$ we get the conditions $a_0>0$, $a_1>0$ and $a_2>0$. For $n=3$ we get the conditions $a_0>0$, $a_1>0$, $a_1a_2-a_0a_3>0$ and $a_3>0$. Note that the third condition is invariant under the replacement of $a_j$ by $a_{n-j}$. When $a_0a_3-a_1a_2>0$, $a_0>0$ and $a_3>0$ then the conditions $a_1>0$ and $a_2>0$ are equivalent to each other. In this way the invariance under reversal of the order of the coefficients becomes manifest. For $n=4$ we get the conditions $a_0>0$, $a_1>0$, $a_1a_2-a_0a_3>0$, $a_1a_2a_3-a_1^2a_4-a_0a_3^2>0$ and $a_4>0$.

Next we look at the issue of loss of stability. If $H$ is the region in matrix space where the Routh-Hurwitz criteria are satisfied, what happens on the boundary of $H$? One possibility is that at least one eigenvalue becomes zero. This is equivalent to the condition $a_n=0$. Let us look at the situation where the boundary is approached while $a_n$ remains positive, in other words the determinant of the matrix remains non-zero. Now $a_0=1$ and so one of the quantities $\Delta_k$ with $1\le k\le n-1$ must become zero. In terms of eigenvalues what happens is that a number of complex conjugate pairs reach the imaginary axis away from zero. The generic case is where it is just one pair. An interesting question is whether and how this kind of event can be detected using the $\Delta_k$ alone. The condition for exactly one pair of roots to reach the imaginary axis is that $\Delta_{n-1}=0$ while the $\Delta_k$ remain positive for $k. In a paper of Liu (J. Math. Anal. Appl. 182, 250) it is shown that the condition for a Hopf bifurcation that the derivative of the real part of the eigenvalues with respect to a parameter is non-zero is equivalent to the condition that the derivative of $\Delta_{n-1}$ with respect to the parameter is non-zero. In a paper with Juliette Hell (Math. Biosci. 282, 162), not knowing the paper of Liu, we proved a result of this kind in the case $n=3$.

### Poincaré, chaos and the limits of predictability

March 5, 2017

In the past I was surprised that there seemed to be no biography of Henri Poincaré. I recently noticed that a biography of him had appeared in 2013. The title is ‘Henri Poincaré. A scientific biography’ and the author is Jeremy Gray. At the moment I have read 390 of the 590 pages. I have learned interesting things from the book but in general I found it rather disappointing. One of the reasons is hinted at by the subtitle ‘A scientific biography’. Compared to what I might have hoped for the book concentrates too much on the science and too little on the man. Perhaps Poincaré kept his private life very much to himself and thus it was not possible to discuss these aspects more but if this is so then I would have found it natural that the book should emphasize this point. I have not noticed anything like that. I also found the discussion of the scientific topics of Poincaré’s work too technical in many places. I would have preferred a presentation of the essential ideas and their significance on a higher level. There are other biographies of great mathematicians which made a better impression on me. I am thinking of the biography of Hilbert by Constance Reid and even of the slim volumes (100 pages each) on Gauss and Klein written in East Germany.

On important discovery of Poincaré was chaos. He discovered it in the context of his work on celestial mechanics and indeed that work was closely connected to his founding the subject of dynamical systems as a new way of approaching ordinary differential equations, emphasizing qualitative and geometric properties in contrast to the combination of complex analysis and algebra which had dominated the subject up to that point. The existence of chaos places limits on predictability and it is remarkable that these do not affect our ability to do science more than they do. For instance it is known that there are chaotic phenomena in the motion of objects belonging to the solar system. This nevertheless does not prevent us from computing the trajectories of the planets and those of space probes sent to the other end of the solar system with high accuracy. These space probes do have control systems which can make small corrections but I nevertheless find it remarkable how much can be computed a priori, although the system as a whole includes chaos.

This issue is part of a bigger question. When we try to obtain a scientific understanding of certain phenomena we are forced to neglect many effects. This is in particular true when setting up a mathematical model. If I model something using ODE then I am, in particular, neglecting spatial effects (which would require partial differential equations) and the fact that often the aim is not to model one particular object but a population of similar objects and I neglect the variation between these objects which I do not have under control and for whose description a stochastic model would be necessary. And of course quantum phenomena are very often neglected. Here I will not try to address these wider issues but I will concentrate on the following more specific question. Suppose I have a system of ODE which is a good description of the real-world situation I want to describe. The evolution of solutions of this system is uniquely determined by initial data. There remains the problem of sensitive dependence on initial data. To be able to make a prediction I would like to know that if I make a small change in the initial data the change in some predicted quantity should be small. What ‘small’ means in practice is fixed by the application. A concrete example is the weather forecast whose essential limits are illustrated mathematically by the Lorenz system, which is one of the icons of chaos. Here the effective limit is a quantitative one: we can get a reasonable weather forecast for a couple of days but not more. More importantly, this time limit is not set by our technology (amount of observational data collected, size of the computer used, sophistication of the numerical programs used) but by the system itself. This time limit will not be relaxed at any time in the future. Thus one way of getting around the effects of chaos is just to restrict the questions we ask by limits on the time scales involved.

Another aspect of this question is that even when we are in a regime where a system of ODE is fully chaotic there will be some aspects of its behaviour which will be predictable. This is why is is possible to talk of ‘chaos theory’- I know too little about this subject to say more about it here. One thing I find intriguing is the question of model reduction. Often it is the case that starting from a system of ODE describing something we can reduce it to an effective model with less variables which still includes essential aspects of the behaviour. If the dimension of the reduced model is one or two then chaos is lost. If there was chaos in the original model how can this be? Has there been some kind of effective averaging? Or have we restricted to a regime (subset of phase space) where chaos is absent? Are the questions we tend to study somehow restricted to chaos-free regions? If the systems being modelled are biological is the prevalence of chaos influenced by the fact that biological systems have evolved? I have seen statements to the effect that biological systems are often ‘on the edge of chaos’, whatever that means.

This post contains many questions and few answers. I just felt the need to bring them up.

### 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$\kappa$B 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$\kappa$B 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$\kappa$B could skew the ratio. I found a video by Baltimore from 2001 about NF$\kappa$B. 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$\kappa$B 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$\kappa$B is activated. There are lots of animated videos on Youtube illustrating this but I prefer a description in words. Normally NF$\kappa$B is found in the cytosol bound to an inhibitor I$\kappa$B. Under certain circumstances a complex of proteins called IKK forms. The last K stands for kinase and IKK phosphorylates I$\kappa$B. This causes I$\kappa$B to be ubiquinated and thus marked for degradation (cf. the discussion of ubiquitin here). When it has been destroyed NF$\kappa$B 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$\kappa$B. What genes does NF$\kappa$B 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$\kappa$B? 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$\kappa$B 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$\kappa$B 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$\kappa$B system. The primary thing which has been observed to oscillate is the concentration of NF$\kappa$B 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$\kappa$B 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$\kappa$B in the nucleus, I$\kappa$B in the cytoplasm and mRNA coding for I$\kappa$B. 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$\kappa$B? 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.

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