Hybrid models of the cell cycle molecular machinery
Piecewise smooth hybrid systems, involving continuous and discrete variables, are suitable models for describing the multiscale regulatory machinery of the biological cells. In hybrid models, the discrete variables can switch on and off some molecular interactions, simulating cell progression through a series of functioning modes. The advancement through the cell cycle is the archetype of such an organized sequence of events. We present an approach, inspired from tropical geometry ideas, allowing to reduce, hybridize and analyse cell cycle models consisting of polynomial or rational ordinary differential equations.
Hybrid systems are widely used in automatic control theory to cope with situations arising when a finite-state machine is coupled to mechanisms that can be modeled by differential equations . It is the case of robots, plant controllers, computer disk drives, automated highway systems, flight control, etc. The general behavior of such systems is to pass from one type of smooth dynamics (mode) described by one set of differential equations to another smooth dynamics (mode) described by another set of differential equations. The command of the modes can be performed by changing one or several discrete variables. The mode change can be accompanied or not by jumps (discontinuities) of the trajectories. Depending on how the discrete variables are changed, there may be several types of hybrid systems: switched systems , multivalued differential automata , piecewise smooth systems . Notice that in the last case, the mode changes when the trajectory attains some smooth manifolds. In these examples, the changes of discrete variables and the evolution of continuous variables are deterministic. The class of hybrid systems can be extended by considering stochastic dynamics of both continuous and discrete variables, leading to piece-wise deterministic processes, switched diffusions or diffusions with jumps [20, 7, 6, 24, 5]. Hybrid, differential, or stochastic Petri nets provide equivalent descriptions of the dynamics and were also used in this context .
The use of hybrid models in systems biology can be justified by the temporal and spatial multi-scaleness of biological processes, and by the need to combine qualitative and quantitative approaches to study dynamics of cell regulatory networks. Furthermore, hybrid modelling offers a good compromise between realistic description of mechanisms of regulation and possibility of testing the model in terms of state reachability and temporal logics [13, 15]. Threshold dynamics of gene regulatory networks [3, 21] or of excitable signaling systems  has been modelled by piecewise-linear and piecewise-affine models. These models have relatively simple structure and can, in certain cases, be identified from data [19, 9].Some methods were proposed for computing the set of reachable states of piecewise affine models .
Among the applications of hybrid modeling, one of the most important is the cell cycle regulation. The machinery of the cell cycle, leading to cell division and proliferation, combines slow growth, spatio-temporal re-organisation of the cell, and rapid changes of regulatory proteins concentrations induced by post-translational modifications. The advancement through the cell cycle is a well defined sequence of stages, separated by checkpoint transitions. This justifies hybrid modelling approaches, such as Tyson’s hybrid model of the mammalian cell cycle . This model is based on a Boolean automaton whose discrete transitions trigger changes of kinetic parameters in a set of ODEs. The model has been used to reproduce flow cytometry data. Instead of building the hybrid model from scratch, another strategy is to identify hybrid models from experimental or artificial time series [17, 18, 2]. The resulting cell cycle hybrid models can be used for hypothesis testing, as they are or as parts of larger, integrated models.
In this paper we develop ideas first introduced in . We discuss how a given model of the cell cycle, based on ODEs, can be hybridized. The hybridization, based on a tropical geometry heuristics, unravels commonalities of cell cycle models. These systems combine quasi-equilibrium states, represented by slow invariant manifolds and excitability, represented by rapid transitions to and from these manifolds. With respect to  we introduce several new concepts and provide rigorous justification of the procedures. Two general hybridization procedures called tropicalizations are introduced in section 2. The tropicalized dynamics is guaranteed to be a good approximation for polynomial or rational systems with well separated terms and that satisfy a condition called permanency. In subsection 2.2 we introduce the tropical equilibration as a method to test permanency. In section 3 we apply these methods to a cell cycle biochemical network model.
2 Tropical geometry and hybridization
2.1 General settings
In chemical kinetics, the reagent concentrations satisfy ordinary differential equations:
Rather generally, the rates are rational functions of the concentrations and read
where , , are multivariate polynomials. Here , , , are nonzero real numbers, and are finite subsets of .
Special case are represented by
where , are Laurent polynomials with positive coefficients, , , are finite subsets of .
In multiscale biochemical systems, the various monomials of the Laurent polynomials have different orders, and at a given time, there is only one or a few dominating terms. Therefore, it could make sense to replace Laurent polynomials with positive real coefficients , by max-plus polynomials .
This heuristic can be used to associate a piecewise-smooth hybrid model to the system of rational ODEs (1), in two different ways.
The two-terms tropicalization was used in  to analyse the dependence of steady states on the model parameters. The complete tropicalization was used for the study of the model dynamics and for the model reduction .
For both tropicalization methods, for each occurrence of the Dom operator, one can introduce a tropical manifold, defined as the subset of where the maximum in Dom is attained by at least two terms. For instance, for , such tropical manifold is made of points, segments connecting these points, and half-lines. The tropical manifolds in such an arrangement decompose the space into sectors, inside which one monomial dominates all the others in the definition of the reagent rates. The study of this arrangement give hints on the possible steady states and attractors, as well as on their bifurcations.
2.2 Justification of the tropicalization and some estimates
In the general case, the tropicalization heuristic is difficult to justify by rigorous estimates, however, this is possible in some cases. We state here some results in this direction. To simplify, let us consider the class of polynomial systems, corresponding to mass action law chemical kinetics:
where are multi-indices, and is a small parameter. So, the right hand side of (6) is a sum of monomials. We suppose that coefficients have different orders in :
where for .
We also suppose that the cone is invariant under dynamics (6) and initial data are positive:
The terms (7) can have different signs, the ones with are production terms, and those with are degradation terms.
From the biochemical point of view, the choice (7) is justified by the fact that biochemical processes have many, well separated timescales. Furthermore, we are interested in biochemical circuits that can function in a stable way. More precisely, we use the permanency concept, borrowed from ecology (the Lotka -Volterra model, see for instance ).
The system (6) is permanent, if there are two constants and , and a function , such that
We assume that and are uniform in (do not depend on) as .
For permanent systems, we can obtain some results justifying the two procedures of tropicalization.
Then the difference satisfies the estimate
where the positive constants are uniform in . If the original system (6) is structurally stable in the domain , then the corresponding tropical systems (4) and (5) are also permanent and there is an orbital topological equivalence between the trajectories and of the corresponding Cauchy problems. The homeomorphism is close to the identity as .
The proof of the estimate (9) follows immediately by the Gronwall lemma. The second assertion follows directly from the definition of structural stability.
Permanency property is not easy to check. In the case of systems (6) we can make a renormalization
and suppose that (8) holds for with uniform in .
We seek for renormalization exponents such that only a few terms dominate all the others, for each -th equation (6) as . Let us denote the number of terms with minimum degree in for -th equation as . Naturally, . After renormalization, we remove all small terms that have smaller orders in as . We can call this procedure tropical removing. The system obtained can be named tropically truncated system.
Let us denote the coefficient of the multi-index . If all then we have the following truncated system
If all , in order to find possible renormalization exponents , it is necessary to resolve a family of linear programming problem. Each problem is defined by a set of pairs such that . We define by
and obtain the system of the following inequalities
The following straightforward lemma gives a necessary condition of permanency of the system (6).
Assume a tropically truncated system is permanent. Then, for each , the -th equation of this system contains at least two terms. The terms should have different signs for coefficients , i.e., one term should be a production one, while another term should be a degradation term.
We call “tropical equilibration”, the condition in Lemma 2.5. This condition means that permanency is acquired only if at least two terms of different signs have the maximal order, for each equation of the system (6). This idea is not new, and can be traced back to Newton.
The tropical equilibration condition can be used to determine the renormalization exponents, by the following algorithm.
Step 1. For each let us choose a pair such that and . The sign of the corresponding terms should be different.
Step 2. We resolve the linear system of algebraic equations
for , together with the inequalities (15).
Notice that although that Step 2 has polynomial complexity, the tropical equilibration problem has an exponential number of choices at Step 1.
Assume that, as a result of this procedure, we obtain the system
One can expect that, in a ”generic” case
We can now state a sufficient condition for permanency. Let us consider the first equation (17) with and let us denote . In this notation, the first equation becomes
Since , one has that is a slow function of time and thus we can suppose that are constants (this step can be rendered rigorous by using the concept of invariant manifold and methods from ). The permanency property can be then checked in an elementary way. All rest points of (19) are roots of . If , is an increasing time function and if , is a decreasing time function. A single root of within is given by
These properties entail the following result:
Equation 19 has the permanence property if and only if
For fixed , in these cases we have
The generic situation described by the conditions (18) lead to trivial “chain-like” relaxation towards a point attractor, provided that we have permanency at each step. This result is the nonlinear analogue of the similar result that monomolecular networks with total separation relax as chains and can only have stable point attractors .
The following theorem describes a less trivial situation, when limit cycles are possible.
Assume holds. If the procedure, described above, leads to the permanency property at each step, where , and the last two equations have a globally attracting hyperbolic rest point or globally attracting hyperbolic limit cycle, then the tropically truncated system is permanent and has an attractor of the same type. Moreover, for sufficiently small the initial system also is permanent for initial data from some appropriate domain and has an analogous attracting hyperbolic rest point (limit cycle) close to the attractor of the truncated system. If the rest point (cycle) is not globally attracting, then we can say nothing on permanency but, for sufficiently small , the initial system still has an analogous attracting hyperbolic rest point (limit cycle) close to the attractor of truncated system and the same topological structure.
Finally, let us note that tropical equilibrations with permanency imply the existence of invariant manifolds. This allows to reduce the number of variables of the model while preserving good accuracy in the description of the dynamics. The following Lemma is useful in this aspect.
Consider the system
where , is a parameter and the function enjoys the following properties. This function lies in an Hölder class
and the corresponding norms are uniformly bounded in , for some open domain :
Assume one of conditions i, ii, iii of Lemma 2.6 holds. We also suppose that are smooth functions of for all such that . Assume that implies .
Proof. The proof is standard, follows from Theorems in , Ch. 9.
3 A paradigmatic cell cycle model and its tropicalization
3.1 Description of the model
We study here the cell cycle model proposed by Tyson . This model mimics the interplay between cyclin and cyclin dependent kinase cdc2 (forming the maturation promoting factor MPF complex) during the progression of the cell cycle. The model demonstrates that this biochemical system can function as an oscillator, or converge to a steady state with large MPF concentration, or behave as an excitable switch. The three regimes can be associated to early embryos rapid division, metaphase arrest of unfertilized eggs, and growth controlled division of somatic cells, respectively. This model takes into account autocatalytic activity of MPF (positive feed-back). It can be described as a nonlinear cycle of biochemical reactions and corresponds to the following set of differential equations:
Here are rate constants, are concentrations of cdc2, p-cdc2 (phosphorylated kinase), cyclin-p:cdc2 complex (active MPF), cyclin-p:cdc2-p complex (inactive MPF), and cyclin, respectively. With respect to the original model we have introduced a small parameter to cope with the order of the rate constants ( in the original model).
The system (26) has the conservation law
where the value (total initial concentration of kinase cdc2) was chosen by convenience.
3.2 Tropical equilibrations and model reduction
Let us apply the tropical equilibration principle. To this aim, we renormalize the variables,
Let us substitute these relations into the system of equations. As a result, we obtain
The system (29) has the conservation law
In order to compute the exponents we use tropical equilibrations together with the conservation law (30). There are variants of tropical equilibrations. To our surprise, there is only one solution for the exponents values. We can show that all possible equilibrations of the variables , and uniquely set the values of two exponents, , .
Let us consider the variants with respect to the equilibrations of the variables and . Denoting by the term in the equation, we have the following situations:
1) In eq. for : , In eq. for : .
2) In eq. for : , In eq. for : .
3) In eq. for : , In eq. for : .
4) In eq. for : , In eq. for : .
In the variant 1 (Case I) the tropical equilibrations do not fix the values of the exponent and we get
In the variants 2,3,4 (Case II) the exponents are uniquely determined from equilibrations and we obtain
However, (32) and the conservation law (30) are incompatible, therefore, the case II can be rejected. In the case I the conservation law takes the form as . Assuming that and (it is reasonable, since it is a ”generic case”), we obtain . Thus, the only possible situation is variant 1 (Case I) and the corresponding set of exponents is:
Let us note that the terms and in the equations for the variables correspond to direct and reverse rates of a phosphorylation/dephosphorylation cycle transforming into and back. Thus, biochemically, (Case I) corresponds to the quasi-equilibrium of this cycle. Furthermore, the equilibration of all the variables leads to the exponents (31). In this case, , , meaning that the tropical equilibrations of the variables , are triple (in each equation, all three terms have the same order).
We finally obtain the following renormalized system
The structure of the system (34) emphasizes the multiple time scales of the model. The fastest variables are in order , then and . The variables and are slow.
This important assumption ensures the existence of an invariant manifold and will be justified, a posteriori.
Then, from the last equation (34) we obtain the relation
which represents the equation of an invariant manifold.
In turn, this relation leads to the following equations for
A second invariant manifold equation is defined by the equation
Remind that is a slow variable. Then, for we have
System (36) represents a two-dimensional reduced model of the initial five-dimensional system. This result shows that tropical equilibrations can be used for model reduction.
The solutions of (36) either tend to the stable equilibrium
or, if this equilibrium is unstable, to a limit cycle.
Based on the general Theorem 2.7 we can assert the following:
Assume (35) holds with . If the shorted system (36) has a stable hyperbolic limit cycle, then, under above conditions, for sufficiently small the five component system (26) also has a stable limit cycle. If the shorted system (36) has a stable hyperbolic equilibrium, then, under above conditions, for sufficiently small the five component system (26) also has a stable hyperbolic equilibrium.
We have studied the system (36) analytically and numerically. The numerical simulations confirm the criteria of cycle existence both for small and for . For small epsilon the cycle has a singular structure. The amplitude of and the cycle period increase in , approximatively, as (the assertion about period is natural since the rate of is ).
Hyperbolicity can be straightforwardly checked for the rest point (39), by computing the eigenvalues of the linearized system. Denote by . When the rest point is hyperbolic and stable we have the following estimate
with some holds. Integrating (38) for over interval gives
uniformly in as . This yields that if and therefore, does not go to zero for large , justifying the estimate (35) needed for the existence of an invariant manifold. The case of a limit cycle is discussed in the next subsection.
3.3 Singular limit cycle and hybrid dynamics
Up to this point, the tropical ideas were used for reducing the dynamics of the model. In this section we show that the tropicalization heuristic is well adapted for decomposing the limit cycle into slow and fast modes, providing a hybrid description of the dynamics.
Let us note that in a hybrid, excitable system, it is possible that not all variables are equilibrated. Also, the system can have more than two different equilibrations and associated invariant manifolds, and jump from one invariant manifold to another during the dynamics. Let us consider that the variables are equilibrated as above, but now, only one among the variables or are equilibrated. We have four situations:
1) In eq. for : , 2) In eq. for : ,
3) In eq. for : , 4) In eq. for : .
Combined with the conservation law condition (30), variants 1 and 2 lead to the same triple tropical equilibration as before (Case I) and exponents (33). We denote the corresponding invariant manifold . The renormalized equations are the same as (36).
Variant 3 can be rejected by the general permanency criterion given by Lemma 2.6. Variant 4 (Case III) satisfies the permanency criterion and leads to
This corresponds to a double equilibration (two equal terms) of the variable , the variable being not equilibrated. We denote the corresponding invariant manifold . The renormalized equations read
We can provide a hybrid description of the cell cycle, by decomposing the periodic orbit into three modes (Fig.1). The first mode is the slowest and has the longest duration. It consists in the dynamics on the slow invariant manifold at low values of , and can be described by the algebraic-differential system (part between and of the orbit in Fig1c)). The exit from the invariant manifold occurs at a critical point (point ). The next slowest mode corresponds to the decrease of (part between and of the orbit in Fig1c)) and can be described by two terms truncated system . Finally, there is a fast mode, corresponding to the fast increase of (part between and of the orbit in Fig1c)) and described by the truncated system . One can notice (Fig1c)) that this hybrid approximation is very accurate for small . At a distance from the tropical manifold, the hybrid orbit coincides with the one generated by the two terms or by the complete tropicalization. However, close to the tropical manifold, the two term and the complete tropicalization are less accurate than the hybrid approximation described above.
Below we state rigorous estimates describing the slow movement on and the fast jump towards . The two terms description of the dynamics on is a direct consequence of the Proposition 2.4 and Lemmas 2.6,2.8.
To simplify notation, we rewrite the system (36) for as
We can obtain such a presentation by a linear variable change.
Let us define the functions
and the points
Let us find some estimates of solutions at . Our goal is to prove that at this point starts to increase sharply. After this, the terms play the main role in the equations (44),(45), and the other terms can be removed while the -component is big.
Let us introduce new variables by
For one obtains
Let us consider for this system the Cauchy problem with initial data
This result can be reinforced. Actually, attains values of the order .
for such that
and such that
We showed that tropical ideas can be usefully employed to reduce and hybridize polynomial or rational dynamical systems occurring in modelling the molecular machinery of the cell cycle. The main idea consists in keeping only the dominant monomial terms in the right hand side of the ordinary differential equations. Depending on the position in phase space, one should keep one, two, or more such terms. The places where two or more monomial terms are equal form the so-called tropical manifolds. The one term approximation is valid far from the tropical manifolds, whereas close to tropical manifolds several dominating terms of opposite signs can equilibrate each other. These “tropical equilibrations” of the dominating terms slow down the dynamics and produce attractive invariant manifolds.
The possible applications of this method are multiple. Generally, the method can be used to obtain simplified models. In the example studied here, we have started with a five variables model, that has been reduced to two variables and hybridized. The modes of the hybrid model have the simple structure of monomial differential or differential-algebraic equations. Two general methods that we called complete and two terms tropicalizations provide description of the modes and of the mode changes. However, these general procedures may lead to inaccurate approximations when the full model does not satisfy permanency globally. In such cases, more thorough analysis is needed. We have shown that the model of embryonic cell cycle has essentially three modes with different timescales, namely slow accumulation of cyclin, rapid activation of MPF and intermediately rapid degradation of cyclin and inactivation of MPF. The fastest mode is described by monomial ODEs, whereas the less fast modes correspond to tropical equilibrations and are described by differential-algebraic equations.
Several improvements and developments are needed in order to apply these methods at a larger scale. The computation of tropical equilibrations suffers from combinatorial explosion. However, for the biochemical network used as working example, the number of solutions seems to be very small compared to the large combinatorics of monomial terms. There is hope, that once formulated in constraint logic programming, the problem of equilibrations could be efficiently computed in practice as a constraint satisfaction problem. Also, effective methods are needed to compute the transitions between modes. The main difficulty here is related to walls (segments of the tropical manifolds) crossing. Near walls, two or more terms are dominant. When these terms are equilibrated, orbits remain close to the walls and are contained in invariant manifolds. The complete or two terms tropicalizations provide general heuristic for mode transitions. These approximations mail fail close to walls. For instance, as we showed in a previous paper , the complete tropicalization predicts sliding modes that evolve on the wall and stay thus close to orbits of the full system. However, these sliding modes can be too long, leaving the wall when the orbits of the full system are already far away. In order to get an accurate description of the behavior near such walls we had to compute invariant manifolds. Although this is generally much simpler than integrating the full set of equations, it could become difficult for tropical equilibrations involving more than two terms. Future work will be dedicated to developing general methods for this problem.
Appendix A Appendix: proofs
Proof of Lemma 3.2 Let us consider the equation
One observes that
Therefore, if , then the solution attains a small - neighborhood of within a bounded time interval