Multiple steady states in a mathematical model for interactions between T cells and macrophages

The aim of this paper is to prove results about the existence and stability of multiple steady states in a system of ordinary differential equations introduced by R. Lev Bar-Or to model the interactions between T cells and macrophages. Previous results showed that for certain values of the parameters these equations have three stationary solutions, two of which are stable. Here it is shown that there are values of the parameters for which the number of stationary solutions is at least seven and the number of stable stationary solutions at least four. This requires approaches different to those used in existing work on this subject. In addition, a rather explicit characterization is obtained of regions of parameter space for which the system has a given number of stationary solutions.


Introduction
Many phenomena in biology and medicine can be modelled by systems of ordinary differential equations. Sometimes models are used which have very large numbers of unknowns and parameters. Often quantitative information about these parameters is obtained from high throughput experimental techniques. This kind of approach plays a central role in the field of systems biology. (For an introduction to systems biology see [3].) A complementary approach is to try to capture important effects with a description in terms of dynamical systems with a few unknowns. At the same time large numbers of parameters may be handled by looking for qualitative features of the solutions which are present for large regions of parameter space. This leads to mathematical problems more accessible to analytical treatment and less dependent on large-scale simulations. Some thoughtful remarks on the relative advantages of these different paths to understanding in the applications of mathematics to biology can be found in [5].
One type of qualitative behaviour which is frequently of interest for biologists is that of multistability, i.e. the existence of several stable steady state solutions for given values of the parameters. This means that a biological system can function in more than one way, depending on its history. In addition to the intrinsic interest of this fact it opens up the possibility of manipulating the biological system in some way which is favourable for practical applications such as therapies for certain diseases.
The subject of what follows is a system of four ordinary differential equations with eighteen parameters introduced in [4] to model certain aspects of the immune system. A mathematically rigorous analysis of some featues of the behaviour of solutions of this system was carried out in [7]. In particular it was proved that for some values of the parameters there is only one stationary solution while for other parameter values there are three stationary solutions, two of which are stable. When there is bistability the two stationary solutions can be identified with two different states of the immune system which are said to be Th1-dominated and Th2-dominated, respectively. They are distinguished by the concentrations of certain substances called cytokines which immune cells use to communicate with each other. In the derivation of the model the cells which produce and react to the cytokines are T cells and macrophages. The unknowns in the ODE system are concentrations of cytokines and the cells do not occur explicitly in the model.
The function of the immune system is to defend the host against harmful influences such as pathogens and toxins which can cause diseases. (For a detailed introduction to immunology see [6] or [8].) Under some circumstances the immune system may malfunction by attacking host tissues, leading to autoimmune diseases. The choice between Th1-and Th2-dominated immune responses may have an important influence on the course of diseases, whether they are due to pathogens or autoimmune in nature. More information on this can be found in [7] and [4]. A survey of work on the development of Th1 and Th2 responses in an individual organism and approaches to understand this phenomenon using mathematical modelling is given in [1].
In what follows it is shown that for suitable values of the parameters the system of [4] has seven stationary solutions, four of which are stable. This means for instance that it is possible to have two steady states which are Th1-dominated with the degree of dominance being different. This property of the system had not previously been observed. As noted in [7], the property of bistability can be obtained if only T cells are taken into account. On the other hand the influence of the macrophages is essential in order to fulfil the assumptions of the theorems in what follows which assert the existence of more than two stable steady states.
The essential qualitative features of the dynamics exhibited in what follows can already be found in a much simpler system which arises in a special case. This model system has only two unknowns and two parameters. The basic system of [4] and the two-dimensional model system are introduced in Sect. 2. Sect. 3 contains statements and proofs of theorems on multistability for the model system. This allows basic ideas to be presented in a relatively simple context. Some of the techniques used extend in a straighforward way to the analysis of the four-dimensional system. In fact they suffice to prove the existence of four stable stationary solutions for an open set of the parameter space which is described rather explicitly. To prove the existence of more steady states of the four-dimensional system a significant refinement of the techniques is necessary. These ideas are presented in Sect. 4. The last section is devoted to further discussion of the results of the paper and possible generalizations.

The basic equations
The basic dynamical system studied in this paper is The d i are positive constants. The function g is given by where θ is a constant. The functions h i are defined by h i = j a ij x j for some constants a ij . The equations were written in this form in [7]. The system was first introduced in [4] in a different notation. The relation between the two notations is explained in [7]. The coefficients in (1) are required to satisfy a number of conditions which will now be listed. The formulation of these conditions has been changed slightly from that used in [7] so as to make certain calculations more efficient. In the formulation used here each a ij is of the form b ij + c ij and the following conditions hold: there is a constant K > 0 such c 2j = Kc 1j for all j.
5. c 11 ≥ 0 and c 13 ≥ 0 The last of these conditions is a slight weakening compared to [4], where strict inequality was required. The analysis done in what follows only requires the assumptions 1. and 3. Conditions 2., 4. and 5. have been included here only to make the relation to the original set-up of [4] clear. The interpretation of these equations is as follows. The quantities x 1 and x 3 are the concentrations of Th1 cytokines produced by T cells and macrophages respectively. The quantities x 2 and x 4 are the corresponding concentrations of Th2 cytokines. Thus z 1 = x 1 + x 3 and z 2 = x 2 + x 4 are the total concentrations of Th1 and Th2 cytokines, respectively. The signs of the b ij reflect the fact that each type of cytokine (Th1, Th2) stimulates cells to produce that type of cytokine and inhibits their production of the other type. The coefficients c ij encode the effect of antigen presentation. The cytokines modify the way in which macrophages present antigen to T cells and this in turn modifies the cytokine production of the T cells. The signs of these coefficients are based on experimental data reviewed in [4]. Setting the coefficients c ij to zero corresponds to neglecting the effect of antigen presentation.
A special case of this system which is important in what follows is obtained by requiring that d i = 1 for all i, θ = 0, the quantities |b ij | with j odd are all equal, the quantities |b ij | with j even are all equal and the quantities c ij with i = 1, 2 are all equal. Then defining A = b 11 , B = b 12 and C = c 11 gives a system which is (21)-(24) of [7]. It implies a closed system for the variables z 1 and z 2 which is (25)-(26) of [7]. The special case A = B of the last system is given by and is a useful model system for gaining intuition about the behaviour of solutions of (1). The constant A is positive while C is non-negative. This system is symmetric under interchange of z 1 and z 2 . Suppose that ( is a stationary solution of (1) with the restrictions on the parameters a ij leading to (3)-(4) and let z * is a stationary solution of (3)-(4) and x * 4 = g(−A(z * 1 − z * 2 )).
Conversely, if (z * 1 , z * 2 ) is a stationary solution of (3)- (4) and (x * 1 , x * 2 , x * 3 , x * 4 ) is defined by (5)-(8) then a stationary solution of (1) is obtained. The fact that the starting solution was stationary assures the consistency conditions that z * Thus proving the existence of stationary solutions of (5)-(8) also leads to a proof of the existence of corresponding stationary solutions of (1).

Existence of stationary solutions
In this section results are proved about the existence of stationary solutions of (3)-(4). The techniques used to do so will be generalized in the next section to prove analogous results about stationary solutions of (1). A picture which helps to provide an intuitive understanding of the results of this section and which played an important role in developing the theorems and their proofs is given as Fig. 1. The assumptions on the coefficients in (1) imply that C ≥ 0 in the system (3)-(4) derived from it. Nevertheless the possibility that C is negative is allowed in this section. The reason for this is that dealing with negative values of C in (3)-(4) yields intuition for the treatment of (1) in the next section where it is allowed that some of the coefficients c ij are negative. For a positive real number with < 1 let Theorem 1 (i) Suppose that η + < 1 is a positive constant and that |C| A ≤ η + . Then if A is sufficiently large there exists an > 0 such that there is a unique stationary solution of (3) and (4) in each of K 1 ( ) and K 2 ( ).
(ii) Suppose in addition that η − > 1 3 and C A ≥ η − . Then for A sufficiently large can be chosen such that there is a unique stationary solution in each K i ( ) with i = 1, 2, 3, 4. Proof Because of the symmetry of the system it is enough to prove the result for K 1 and K 3 . Define a mapping by Alternatively we can write φ(z 1 , z 2 ) = (φ 1 (z 1 , z 2 ), φ 2 (z 1 , z 2 )). Stationary solutions of (3) and (4) are in one to one correspondence with fixed points of the mapping φ. To prove the theorem it is enough to show that for A sufficiently large and a suitable choice of the mapping φ leaves each K i invariant and that the restriction of φ to each of these sets is a contraction. Suppose that ≤ 1 4 (1 − η + ). Let A 0 be a positive constant and assume that A ≥ A 0 . The hypotheses of the theorem imply that Similarly Next the arguments of g occurring in the second component will be estimated. and Hence If A 0 is large enough then can be chosen so that the inequality e −2(1−η+)A0 ≤ 2 is satisfied. Under these circumstances it follows that φ 1 (z 1 , z 2 ), which is evidently positive, is less than . In the same way g( It follows that with the restrictions on A 0 and already assumed φ 2 , which is evidently less than two, is no less than 2 − . Hence It follows by the mean value theorem that if x and y both have the same sign and modulus greater than x 0 then Hence if (z 1 , z 2 ) and (w 1 , w 2 ) both belong to K 1 then for i = 1, 2. The function xe −2(1−η+)x has its unique maximum on the interval [0, ∞) at A similar analysis can be carried out in the case of K 3 when the assumptions of part (ii) of the theorem hold. For (z 1 , z 2 ) ∈ K 3 .
If it is assumed that e −A0 < 2 and e −(3η−−1)A0 < 2 these inequalities imply that K 3 is invariant. That the restriction of φ to K 3 is a contraction follows as in the case of K 1 under the assumptions that A 0 ≥ 1 3η−−1 This completes the proof of the theorem.
The symmetry of the system implies that the diagonal z 1 = z 2 is an invariant submanifold. Stationary solutions on the diagonal are in one to one correspondence with points where d dt (z 1 + z 2 ) = 0. When Thus stationary solutions on the diagonal are in one to one correspondence with zeroes of the function F (x) = 1 + 2g(Cx) − x on the interval (0, 4). The facts that F (0) = 2 and F (4) < −1 show that F has a zero. The function x − 1 is monotone increasing while for C ≤ 0 the function g(Cx) is monotone nonincreasing. Thus in that case the zero of F is unique. For C > 0 on the other hand, when the graph of 2g(Cx) crosses that of x − 1 its derivative must be no greater than one. Moreover in this case the second derivative of 2g(Cx) is negative and so the derivative of 2g(Cx) is less than one for all greater values of x. Thus its graph cannot cross that of x − 1 again. It follows that the zero of F is unique in all cases and that there is always precisely one stationary solution on the diagonal. Hence in the cases where Theorem 1 shows the existence of four stationary solutions there are at least five in total. Note that if C is large and positive then the stationary point is close to ( 3 2 , 3 2 ). Theorem 2 Under the hypotheses of Theorem 1 the system (3)-(4) has at least seven stationary solutions. Proof By Theorem 1 and the discussion following it there exist at least five stationary solutions. Using the symmetry of the system it is enough to show that there exists at least one other stationary solution. Consider the nullclines of the system (3)-(4), i.e. the zero sets of the right hand sides of the equations. On the first nullcline if z 1 is fixed the monotonicity of g implies that there is at most one possible value for z 2 . Combining this with the implicit function theorem shows that the nullcline is the graph of a function γ 2 of z 1 defined on some open subset. On any closed subinterval of this subset the value of z 2 is bounded. It follows that in fact γ 2 is defined on the whole interval (0, 2). In a similar way it can be shown that the second nullcline can be written in the form z 1 = γ 1 (z 2 ) for a function γ 1 defined on (0, 2). The part of the curve z 1 = γ 1 (z 2 ) where z 1 ≤ 1+ is contained in the strip defined by 2 − 2e −2A0 ≤ z 2 < 2. This can be seen as follows. From estimates very similar to those derived in the proof of Theorem 1 it can be concluded that for any fixed value of z 1 no greater that 1 + the mapping z 2 → φ 2 (z 1 , z 2 ) defined on the interval [2 − , 2] is a contraction. The unique fixed point of this mapping lies on the curve z 1 = γ 1 (z 2 ) and satisfies the estimate claimed. Let z * 1,i and z * 2,i be the z 1 and z 2 coordinates of the stationary solution in K i , respectively. To prove the theorem it is enough to show that in the interval (z * 1,1 , z * 1,3 ) there is a value of z 1 for which γ 2 is less than 2 − 2e −2A0 and a value of z 1 for which γ 2 is greater than 2. For this it is helpful to compute the derivative of γ 2 .
The second expression in square brackets is small and, in particular, less than one half. Taking it onto the left hand side shows that dγ2 dz1 is negative and that It can be concluded that in K 1 and K 3 . This guarantees the existence of values of z 1 with the desired properties and completes the proof. When the linearization of the system about a stationary point in one of the K i , 1 ≤ i ≤ 4, is computed the contributions of the nonlinear terms are extremely small. As a consequence the linearization about a point of this type is very close to minus the identity. It follows that the stationary point is a hyperbolic sink.

Stationary solutions of the full system
In this section some of the techniques introduced in Sect. 3 are applied to obtain information about stationary solutions of (1) which covers an open set of parameters for this system whose definition is based on explicit inequalities. The system of equations defining stationary solutions can be simplified by replacing the variables x i by d i x i . This means that when proving theorems it can be assumed without loss of generality that d i = 1 for all i. Corresponding results for general values of d i can then be obtained by transforming the stationary solutions back to the original variables. First Theorem 1 will be generalized. Note that if a stationary solution of (3)-(4) is close to (0, 2), (2, 0), (1, 2) or (2, 1) then the corresponding stationary solution of (1) defined in Sect. 2 is close to (0, 1, 0, 1)), (1, 0, 1, 0), (1, 1, 0, 1) or (1, 1, 1, 0) respectively. For a positive real number with < 1 definẽ Theorem 3 (i) Suppose that d i = 1 for each i, that η + < 1 and ∆ are positive constants, that |c ij | ≤ η + |b ij | for all i and j and that max i,j {|b ij |} ≤ ∆ min i,j {|b ij |}. Then if min i,j {|b ij |} is sufficiently large there exists an > 0 such that there is a unique stationary solution of (1) in each ofK 1 ( ) andK 2 ( ). (ii) Suppose in addition that ∆−1 is sufficiently small, η − > 1 3 , c ij ≥ η − max k,l {|b kl |} for i = 1, 2 and all j. Then for min i,j {|b ij |} sufficiently large there exists an such that there is a unique stationary solution in eachK i ( ) with i = 1, 2, 3, 4.
Proof The proof uses the same strategy as that of Theorem 1. Define a mapping φ : Stationary solutions of (1) are in one to one correspondence with fixed points of φ. To prove the theorem it suffices to show that each setK i is invariant under φ and that the restriction ofφ to eachK i is a contraction. In comparison with the proof of Theorem 1 the estimates here give less explicit information about the conditions in the statement of the theorem. Let A 0 be a positive constant and assume that |b ij | ≥ A 0 for all i and j. Consider first the setK 1 . Then and this is the case if A 0 is sufficiently large and is chosen appropriately. It can be shown that the restriction ofφ toK 1 is a contraction using a straightforward generalization of the argument used in the proof of Theorem 1. Up to a numerical factor the contraction constant is and this can be made as small as desired by choosing A 0 large enough. The same arguments can be used to prove thatK 2 is invariant underφ and that the restriction ofφ toK 2 is a contraction. Next part (ii) of the theorem will be proved. This will be presented for the case ofK 3 -the case ofK 4 can be treated in the same way. Note that if ξ = 1 − ∆ −1 then assuming that ∆ − 1 is sufficiently small is equivalent to assuming that ξ is sufficiently small. Moreover for all i, j, k, l the quantity ||b ij | − |b kl || is bounded by ξ max i,j |b ij |. In this case This shows that for ξ and sufficiently small the setK 3 is invariant underφ. That the restriction toK 3 is a contraction can be shown as in the case ofK 1 .
Note that for the parameter values which reduce (1) to (3)-(4) the hypotheses of Theorem 3 reduce to those of Theorem 1.
The argument used to obtain the stationary solution on the diagonal for the system (3)-(4) does not generalize to the system (1) since the latter has no symmetry. Instead the observation concerning the approximate position of that solution for A large made in the last section will be used. Definẽ The direct analogue of the arguments used to treatK i for 1 ≤ i ≤ 4 does not seem to work forK 5 which is apparently not invariant underφ. To overcome thisφ will be replaced by another mapping ψ. Theorem 4 Assume that the hypotheses of Theorem 3 hold. Then if ∆ − 1 is sufficiently small and min i,j |b ij | is sufficiently large there exists an > 0 such that there is a stationary solution of (1) inK 5 . Proof Define new variables by y 1 = x 1 − 1, y 2 = x 2 − 1, y 3 = x 3 − 1 2 and The first two terms in each of these expressions are obstructions to showing that φ mapsK 5 into itself. To go further note that for any positive real numbers (x, w) This suggests rewriting the equation x i = g(h i ) for i = 3, 4 as where For i = 1, 2 the equation is rewritten as y i = s i where s 1 = g(h 1 ) − 1. These equations are schematically of the form N y = s for some matrix N . Conditions will now be given which ensure that the matrix N is invertible. For this it is useful to write it as a matrix I 0 N 21 N 22 of two-by-two blocks. If the matrix N 22 is invertible then N is also invertible and the inverse of N is given by Now consider the matrix N 22 . Its determinant is Let β = min i,j |b ij |. Then b −1 33 ≤ β −1 . It follows in particular that if β ≥ 2 then 1 − b −1 33 ≥ 1 2 . In addition b33+b43 b33 ≤ ξ. It follows that the square bracket in (35) is bounded above by −2(1 − ξ) and from below by −4(1 + ξ). In can be concluded that the following inequality holds In particular, the determinant does not vanish and so N is invertible. Moreover the modulus of the inverse of the determinant can be bounded by a constant multiple of β −1 . It follows that the norm of N −1 22 , can be bounded by a constant and hence that the norm of N −1 can be bounded by a constant multiple of β. On the translate ofB 5 by −1, −1, − 1 2 , − 1 2 define a mapping by ψ(y) = N −1 s(y).
Choose ξ < 1 2 η − and < 1 10 η − . Then the first bracket on the right hand side of this equation is greater than 1 2 η − . This in turn gives exponential bounds for s 1 and s 2 . It remains to estimate s 3 and s 4 . Only the argument for s 2 will be given since the corresponding argument for s 4 is very similar.
Using this and (33) gives the estimate By choosing ξ and β 2 sufficiently small it can be ensured that s 3 and s 4 are bounded by an arbitarily small multiple of . Now can be chosen to satisfy this condition and to ensure at the same time that s 1 and s 2 are bounded by an arbitrarily small multiple of . In this way it can be ensured that ψ leaves B 5 invariant. The conclusion of the theorem follows by applying the Brouwer fixed point theorem to ψ. Remark This theorem also holds, with the same proof, if the hypothesis on η − is weakened to η − > 0. It has proved possible to extend the method used to obtain Theorem 4 to obtain the existence of analogues of the remaining stationary solutions of (3)- (4). Unfortunately this has only been achieved for a quite restricted choice of parameters. To understand where the condition for the parameters comes from, note that if then g(h 1 (ζ 1 , 1, 0, 1)) = g(h 2 (1, ζ 2 , 1, 0)) = 1 2 (42) The cases treated in the next theorem are those where ζ 1 or ζ 2 is close to 1 2 . LetK Theorem 5 Assume that the hypotheses of Theorem 4 hold and that If ξ is small enough and β is sufficiently large then can be chosen so that there is a stationary solution of (1) inK 6 ( ). An analogous statement holds if the inequality (43) is replaced by andK 6 ( ) is replaced byK 7 ( ) in the conclusion.
Proof The strategy of the proof follows that of Theorem 4. It is enough to treat the case ofK 6 ( ) since that ofK 7 ( ) is very similar. This time define new variables by y 1 = x 1 − 1 2 , y 2 = x 2 − 1, y 3 = x 3 and y 4 = x 4 − 1.
(b 11 + c 11 )x 1 + (b 12 + c 12 )x 2 + (b 13 + c 13 )x 3 + (b 14 + c 14 )x 4 = (b 11 + c 11 )y 1 + (b 12 + c 12 )y 2 + (b 13 + c 13 )y 3 + (b 14 + c 14 )y 4 In this case the equation x 1 = g(h i ) will be rewritten as (47) For i = 2, 3, 4 the equation x i = g(h i ) is written as y 2 = g(h 2 ) − 1, y 3 = g(y 3 ) and y 4 = g(h 4 ) − 1. So as in the proof of Theorem 4 the condition for stationary solutions has been rewritten in the form N y = s for a matrix N . In this case N is of the form N 11 N 12 0 I where N 11 is scalar with the only entry 1 − 1 2 a 11 , N 12 is a row matrix with entries − 1 2 a 1i , i = 2, 3, 4, and the identity matrix is three by three. 1 − 1 2 a 11 ≤ 1 − 1 2 β. If β > 3 the matrix is invertible. The inverse of N 11 can be bounded by a constant multiple of β −1 . The inverse matrix N −1 is easily computed. Its norm can be bounded by a constant multiple of β. On the translate ofB 6 by − 1 2 , −1, 0, −1 define a mapping by ψ(y) = N −1 s(y). It will now be shown that ψ maps its domain into itself. The estimates suffice to take care of the last three components s i . It remains to treat the first component. Once this has been done the estimate for N −1 which is already available completes the proof. The quantity s 1 can be estimated using the inequality (32).
It follows that This estimate leads to the desired conclusion as in the proof of Theorem 4.
For the special choices of the parameters leading to the system (3)-(4) the quantities ζ 1 and ζ 2 reduce to In this case, as C A varies from one to 1 3 , the ζ i vary from zero to one. Each of the extra conditions (43) and (44) on the parameters reduces to C A = 3 5 . When Theorem 5 applies this gives some information about the position of the stationary solutions whose existence is guaranteed by Theorem 2.

Conclusions and outlook
The dynamical system introduced in [4] to describe the interactions between T cells and macrophages has been analysed with respect to its stationary solutions. For certain open sets of the parameter space it was shown that there are at least seven stationary solutions of which at least four are stable. The positions of these stationary solutions were described approximately. The same features are already found in the model system (3)-(4) with the two parameters A and C.
There it is only necessary to assume that 1 3 < C A < 1 and that A is sufficently large to get these conclusions. Under weaker assumptions, which reduce to |C| A < 1 and A sufficiently large for the model system it was shown that there are at least three stationary solutions of which at least two are stable. By contrast it was shown in [7] that under a suitable smallness assumption on the parameters there is a unique stationary solution and that all solutions converge to it at late times. For the model system a sufficient condition implying the latter behaviour is 2A + |C| < 1.
As already remarked in [7], for C < A the system (3)-(4) is competitive and hence any solution converges to a stationary solution as a consequence of the results of [2]. In particular there can be no periodic solutions. This conclusion depends on the fact that the system is two-dimensional and so cannot be drawn for the system (1), even in the case of parameter values for which it is competitive. The results of this paper do not rule out the possibility that there are values of the parameters for which the model system has periodic solutions. They also do not prove that there cannot be more than seven stationary solutions or more than four stable stationary solutions for some values of the parameters. Simulations have given no indications that either of these phenomena (periodic solutions or extra stationary solutions) actually occur.
From the biological point of view one of the most interesting new findings of this paper is the coexistence of stationary solutions such as those close to (0, 1, 0, 1) and (1, 1, 1, 0) with the same values of the parameters. Certain small changes of the parameters can move the system from a situation where both of these exist to a situation where one of them has disappeared. On a heuristic level this gives a scenario where a small external influence on the system leads to a large change in its properties. Note that these two stationary solutions have very different biological properties. In the second case the Th1 and Th2 cells are secreting a comparable large amount of cytokines. The total cytokine concentration is 50% greater in the second case than in the first. The model of [4] is crude in may ways but their are many models for biological systems which share qualitative features with this one and the kind of switching phenomenon which has been observed here may occur much more widely. Thus it is good to have mathematical tools which allow it to be studied.
The Th1/Th2 picture of autoimmune diseases is no longer up to date. Newer models involve players other than Th1 and Th2 cells. In particular Th17 cells have come to play a very important role [9]. Mathematical insights obtained for one model can be useful in understanding other models, even those which remain to be invented. Thus the fact that the biological content of the model of [4] is not close to the explanations of the behaviour of the immune system which are presently most popular does not prevent it being worthwhile to study its mathematical properties in detail.
The mathematical description of the immune system in the model of [4] is at the level of the interaction of populations of cells and what happens in a given cell is handled as a simple black box. Another very interesting task is to model the processes within cells which determine whether they differentiate into types such as Th1 and Th2. These usually involve the behaviour of transcription factors. Some models of this kind are studied in [11] and [1]. Models of the population of cells and models of the molecular machinery within cells describe systems whose intrinsic nature are very different. Nevertheless they may be related mathematically and exploiting this kind of relation is a particular strength of mathematics.