Gain and Loss of Function mutations in biological chemical reaction networks: a mathematical model with application to colorectal cancer cells
GGain and Loss of Function mutations in biologicalchemical reaction networks: a mathematical modelwith application to colorectal cancer cells
Sara Sommariva · Giacomo Caviglia · Michele PianaAbstract
This paper studies a system of Ordinary Differential Equationsmodeling a chemical reaction network and derives from it a simulation toolmimicking Loss of Function and Gain of Function mutations found in cancercells. More specifically, from a theoretical perspective, our approach focuseson the determination of moiety conservation laws for the system and theirrelation with the corresponding stoichiometric surfaces. Then we show thatLoss of Function mutations can be implemented in the model via modificationof the initial conditions in the system, while Gain of Function mutations canbe implemented by eliminating specific reactions. Finally, the model is utilizedto examine in detail the G1-S phase of a colorectal cancer cell.
Keywords
Reaction kinetics - Synthetic cell biology - Loss of function mu-tations - Gain of function mutations - Colorectal cancer cells - G1-S transitionpoint
Signalling networks are Chemical Reaction Networks (CRNs) consisting of aninterconnected set of pathways, modeling the flow of chemical reactions initi-
S. SommarivaDipartimento di Matematica, Universit´a di Genova, via Dodecaneso 35 16146 Genova, ItalyTel.: +39-010-3536644E-mail: [email protected]. CavigliaDipartimento di Matematica, Universit´a di Genova, via Dodecaneso 35 16146 Genova, ItalyTel.: +39-010-3536830E-mail: [email protected]. PianaDipartimento di Matematica, Universit´a di Genova, and CNR - SPIN GENOVA, via Dode-caneso 35 16146 Genova, ItalyTel.: +39-010-3536939E-mail: [email protected] a r X i v : . [ q - b i o . M N ] A p r Sara Sommariva et al. ated by information sensed from the environment through families of receptorligands (Jordan et al., 2000; Sever and Brugge, 2015; Tyson and Novak, 2008).The reactions in the network follow the process of information transfer insidecytosol down to the description of the activity of a few related transcriptionfactors (Karin and Smeal, 1992; Kohn et al., 2006). Similarly to what happensfor many networks of biological interest, a signaling network may comprisehundreds of reacting chemical species and reactions.Any CRN can be mapped into a system of ordinary differential equations(ODEs) by following standard procedures (Feinberg, 1987; Yu and Craciun,2018; Chellaboina et al., 2009) thus allowing for simulations of the kineticsof the signaling process in the biologic case (see, e.g., Anderson et al. (2019);Roy and Finley (2017)). In principle, the resulting mathematical model is ca-pable of describing healthy physiologic states, and can be adapted to simulateindividual pathological conditions associated to mutations. Since most cancerdiseases result from an accumulation of a set of mutations (Stratton et al.,2009; Levine et al., 2019), the numerical solution of this mathematical modelrepresents a convenient computational tool for the simulation of the mech-anisms giving rise to a tumor. Further, this kind of models can typically betuned in order to mimic the effects of targeted therapies (Facchetti et al., 2012;Levine, 2019) and drug resistance mechanisms (Eduati et al., 2017).The present paper first realizes an analysis of the system of ODEs asso-ciated to a CRN, characterized by a level of generality sufficient to providean efficient simulation tool. From a formal viewpoint, we apply mathematicalmethods devised for the investigation of deterministic homogeneous chemicalreaction systems based on mass-action kinetics. However, in our approach acrucial role is played by moiety conservation laws (CLs), which are essentialin the determination and interpretation of results (De Martino et al., 2014;Shinar et al., 2009). Further, a geometric classification of equilibrium statesin terms of stoichiometric surfaces is discussed, together with their stabilityproperties.From a more operating viewpoint, we define projection operators, whichmap the system of ODE for physiologic conditions into the mutated systemwhich models Loss of Function (LoF) and Gain of Function (GoF) mutations(Griffiths et al., 2005; Li et al., 2019). The model is built and made operativein such a way that it can be modified almost straightforwardly by additionor elimination of chemical reactions or chemical species, change in the valuesof the rate constants in the formulation of mass-action laws, or change of theinitial conditions. As a consequence, the parameters values that have beenoriginally considered as fixed can be customized to fit the tumor data of aspecific patient.As an application, we examine in detail the kinetics of a recently proposednetwork which simulates how colorectal cancer (CRC) cells process informa-tion from external growth factors and the related answer (Tortolina et al.,2015). The network is focused on the G1-S transition point. In this transitiona newborn cell in the G1 phase of its cycle must pass the control of this check- itle Suppressed Due to Excessive Length 3 point before starting the S phase of synthesis of new DNA (Tyson and Novak,2008); indeed, this is the first necessary step towards proliferation.The structure of the paper is as follows. Section 2 provides the mathemati-cal background for the modeling of CRNs for cell signaling. Section 3 containsthe mathematical model describing LoF and GoF mutation processes. Section4 is devoted to the application of the model to CRC cells. Our conclusions areoffered in Section 5.
We consider a CRN consisting of r reactions, denoted as R j , j = 1 , . . . , r , thatinvolve n well-mixed reacting species, denoted as A i , i = 1 , . . . , n . The networkis modeled as a dynamical system with state vector x = ( x , . . . , x n ) T ∈ R n + ,where the upper T denotes transposition, R + is the set of non-negative realnumber, and the generic component x i is the molar concentration (nM) ofthe species A i . According to this chemical interpretation, x is also called con-centration vector (Yu and Craciun, 2018). We assume that the law of massaction holds: when two or more reactants are involved in a reaction, the reac-tion rate is proportional to the product of their concentrations. The resultingpolynomial system of ODEs for the state variables is written as˙ x = S v ( x , k ) , (1)where the superposed dot denotes the time derivative; k = ( k , . . . , k r ) T ∈ R r + stands for the set of rate constants; S is the R n × r constant stoichiometricmatrix, v ( x , k ) ∈ R r + is the vector of reaction fluxes. Here, system (1) accountsfor internal reaction and boundary fluxes (Kschischo, 2010; Schilling et al.,2000) and hence is open (Feinberg, 1987). The matrix element S ij is the netnumber of molecules of the species A i that are produced whenever the reaction R j occurs. Thus the columns of S have been referred to as reaction vectors.We are now interested in investigating the general properties of the solu-tions of the system (1), and, in particular, in computing the correspondingasymptotically stable states. Our study is based on the analysis of the conser-vation laws and the stoichiometric compatibility classes of the system revisedin the next two subsections.2.1 Conservation laws and elemental species Definition 1
Let x ( t ) be a solution of the system of ODEs (1). A constantvector γ ∈ N n \ { } is said to be a semi–positive conservation vector if thereexists c ∈ R + such that γ T x ( t ) = c ∀ t . (2)Moreover, the relation (2) is called a semi–positive conservation law , or equiv-alently moiety conservation law . Sara Sommariva et al.
Since concentrations are expressed in nM, the biochemical interpretationof the conservation law is that the total number of molecules involved in thespecies combination γ T x remains constant during the evolution in time of thenetwork (Shinar et al., 2009). The constant total amount of available moleculesis fixed by γ and the initial state x through c = γ T x . Moreover, the con-centrations of the species involved in the conservation law are bounded fromabove by the constant c (De Martino et al., 2014; Haraldsd´ottir and Fleming,2016; Shinar et al., 2009; Schuster and H¨ofer, 1991). The following propositionrelates the conservation vectors to the properties of the stoichiometric matrix. Proposition 1 If γ ∈ ker( S T ) ∩ N n \ { } , then γ is a conservation vector.Proof The thesis follows by observing that for any solution x ( t ) of the systemof ODEs (1) ddt ( γ T x ) = γ T ˙ x = γ T Sv ( x , k ) = 0 , where the last term is equal to zero as γ ∈ ker( S T ).We concentrate on conservation vectors γ ∈ ker( S T ). Proposition 1 impliesthat the set of all possible semi–positive conservation vectors defines a convexcone whose independent generators can be computed e.g. through the algo-rithm proposed by Schuster and H¨ofer (1991). In the following we will denotewith (cid:8) γ , . . . , γ p (cid:9) a set of such generators and define the matrix N ∈ R p × n as N = γ T ... γ Tp . (3)Henceforth we will assume p = n − rank( S ) and thus the set (cid:8) γ , . . . , γ p (cid:9) defines a basis for ker( S T )Denote by x the initial point of a trajectory. Since NS = 0, it follows that Nx ( t ) = Nx =: c (4)is the constant vector in R p whose components are the constants involvedin the CLs. Thus the matrix N determines p linear, independent CLs; ac-cordingly, the representative point x ( t ) is constrained to move on the affinesubspace of R n which is determined by equation (4), and hence is identified by N and the initial data. Moreover, the linear affine subspace is the intersectionof p hyperplanes.We conclude with a few remarks that will play a fundamental role in theanalysis of the solutions of system (1). In general, there are chemical speciesthat are not involved explicitly in the expressions of the CLs while the remain-ing species may belong to more than one CL. itle Suppressed Due to Excessive Length 5 Definition 2
We say that a CRN is elemented if: (i) it admits a set of gener-ators (cid:8) γ , . . . , γ p (cid:9) such that the matrix N contains at least one minor equalto the identity matrix of order p , say I p ; (ii) each chemical species is involvedin at least one CL, i.e. for each i = 1 , . . . , n there exists k ∈ { , . . . , p } suchthat γ ki (cid:54) = 0. If only condition (i) is fulfilled, the CRN is weakly elemented .Borrowing the terminology of Shinar et al. (2009), the species associatedwith the minor equal to the identity matrix will be called elemental or basicspecies . Remark 1
Definition 2 implies that each elemental species belongs to one, andonly one conservation law. The idea is that elemental species consist of proteinsin free form, which bind to other species involved in the same conservation law,in order to form the derived compounds or secondary species.In general the set of elemental species of a network is not unique as itdepends on the choice of the basis of the conservation vectors and, fixed abasis, multiple minors equal to the identity matrix may exist.
Remark 2
Given an elemented CRN, up to a change in the order of the com-ponents of x , the matrix N may be decomposed as N = [ I p , N ] (5)with N ∈ R p × ( n − p ) . Similarly, the concentration vector x may be decomposedas x = (cid:18) x x (cid:19) (6)with x ∈ R p , and x ∈ R n − p . In particular, x ∈ R p is formed by theelemental variables. Thus equation (4) may be rewritten in the equivalentform x = c − N x (7)Henceforth we assume the CRN to be weakly elemented, with elementalvariables x ... x p , so that the matrix N is described by the block decomposition(5).According to (7), the set of p conservation equations is solved straight-forwardly with respect to the elemental variables, thus yielding a parametricdescription of the affine space defined in (4). Furthermore, substitution of theexpression (7) of x into the system (1) provides a reduced formulation of theoriginal system of ODEs. Definition 3
Consider an elemented CRN. A concentration vector x ∈ R n is an ideal state for the network iff only the elemental species have non-zeroconcentration, i.e., referring to equation (6), x = 0. Sara Sommariva et al. x (0) = x , consider the corresponding solution x ( t )of the system (1), defined at least in the time domain [0 , T ]. Integration in timeof both sides of (1) leads to x ( t ) − x = S (cid:90) t v ( x ( τ ) , k ) dτ , t ∈ [0 , T ] . (8)This shows that x ( t ) − x is a linear combination of the reaction vectors, withtime-dependent coefficients. Therefore x ( t ) − x belongs to a vector space ofdimension equal to rank( S ) defined by the image of the stoichiometric matrix(Feinberg, 1987, 1995; Yu and Craciun, 2018). Accordingly, we provide thefollowing definition of stoichiometric compatibility class (SCC). Definition 4
Given a value x ∈ R n + of the state variable for system (1), wedefine a stoichiometric compatibility class (SCC) of x the set SC ( x ) = ( x + span( S )) ∩ R n + . (9) Proposition 2
Let (cid:8) γ , . . . , γ p (cid:9) be a set of generators of the the convex conedefined by the semi–positive conservation laws and let N be the matrix definedas in equation (5). If p = n − rank( S ) then for all x ∈ R n + SC ( x ) = (cid:0) x + ker( S T ) ⊥ (cid:1) ∩ R n + (10)= (cid:8) y ∈ R n + s.t. Ny = Nx (cid:9) . (11) Proof
Equation (10) simply follows from Definition 4 by observing that span( S ) =(ker( S T )) ⊥ . Equation (11) follows from the fact that, since p = n − rank( S ),the set (cid:8) γ , . . . , γ p (cid:9) is a basis for ker( S T ). Definition 5
A CRN is said to satisfy the global stability condition if forevery stoichiometric compatibility class there exists a unique globally asymp-totically stable state x e .This means in particular that every trajectory with initial point on thegiven SCC tends asymptotically to the steady state x e , which is also an equi-librium point. Details about asymptotic stability properties of CRNs can befound, e.g., in (Chellaboina et al., 2009; Feinberg, 1987; Yu and Craciun, 2018). A mutation consists essentially in a permanent alteration in the nucleotidesequence of the genome of a cell. Mutations play a fundamental role in cancerevolution (Stratton et al., 2009; Weinstein et al., 2013). Here we are con-cerned with effects induced by mutations on species concentrations in theCRN. Specifically, we consider LoF and GoF mutations that are commonly itle Suppressed Due to Excessive Length 7 observed in cancer cells (Hochman et al., 2017; Lemieux et al., 2015; Levine,2019; Levine et al., 2019).LoF mutations result in reduction or abolishment of a protein function,which is simulated by a restriction on the value of the related density. Thedegree to which the function is lost can vary; for null mutations the functionis completely lost and the concentration of the related molecules is supposedto vanish; for leaky mutations some function is conserved and the value ofconcentration is appropriately reduced (Griffiths et al., 2005; Li et al., 2019).In this study, we deal with null LoF mutations by referring directly to theconcentrations of the mutated molecular species, and by assuming that theyare set equal to zero.GoF mutations are responsible for an enhanced activity of a specific pro-tein, so that its effects become stronger (Griffiths et al., 2005; Li et al., 2019).In our framework a GoF is implemented by excluding from our CRN the re-actions involved in the deactivation of the considered protein; this is achieved,by setting to zero the corresponding reaction rates.3.1 Loss of FunctionConsider a CRN described by system (1) and consider a state x of the system.A LoF mutation of the elemental species A j results in a projection of the state x in a novel state where the concentrations of the j − th elemental species andof all its compounds are zero. Equivalently, also the total concentration c j available in the j − th conservation law is zero. This is modeled by applyingthe following operator to the state x . Definition 6
Consider an elemented CRN and let A j be the j − th elementalspecies of the network. A LoF mutation of A j results in the operator P L j : R n → R n that projects x into a novel state P L j ( x ) = (cid:2) ˜ x T , ˜ x T (cid:3) T such that˜ x ,i = (cid:40) γ ji (cid:54) = 0 x ,i otherwise , i = p + 1 , . . . , n and ˜ x = ( (cid:101) c − N ˜ x )where (cid:101) c is obtained by setting to 0 the j − th element of the vector c = Nx . Remark 3 If x is an ideal state for the CRN then P Lj ( x ) is obtained by settingto zero the concentration x j of the j − th elemental species. Remark 4 P L j transforms SCCs into SCCs. Indeed, if the concentration vec-tors x and y belong to the same SCC, i. e. Nx = Ny = c , Sara Sommariva et al. then P L j ( x ) and P L j ( y ) still belong to the same SCC. More in details, N P L j ( x ) = N P L j ( y ) = ˜ c , Where ˜ c is equal to c except for the j − th element that is zero.Therefore the following theorem holds. Theorem 1
Consider an elemented CRN described by the system of ODEs(1), and any states x , y such that SC ( x ) = SC ( y ) . Then SC ( P L j ( x )) = SC ( P L j ( y )) .If in addition the CRN satisfies the global stability condition then the trajecto-ries starting from P L j ( x ) and P L j ( y ) lead to the same globally asymptoticallystable state.Remark 5 This holds in particular if x and y are replaced by the initial state x and the corresponding steady state x e , respectively. Thus the trajectoriesstarting from P L j ( x ) and P L j ( x e ) lead to the same mutated steady state.3.2 Gain of FunctionConsider a CRN described by system (1). In particular, let S be the stoi-chiometric matrix of the system. A mutation resulting in the GoF of a givenprotein is implemented by removing from the network the reactions involvedin the deactivation of such a protein. From a mathematical viewpoint this canbe achieved by setting to zero the values of the corresponding rate constants,or equivalently by setting to zero the corresponding columns of S . Definition 7
Consider a CRN and let S be the corresponding stoichiometricmatrix. Given a set of reactions identified by the indices H ⊆ { , . . . , r } a GoFmutation results in the operator G H : R n × r → R n × r that projects S into anovel stoichiometric matrix (cid:101) S = G H ( S ) such that˜ S i,h = (cid:40) h ∈ HS i,h otherwise , i = 1 , . . . , n , Remark 6 G H ( S ) defines a new CRN where the chemical reactions in H havebeen removed, while the set of chemical species, and thus the state space R n ,are kept fixed. Theorem 2
If the set of reactions H ⊆ { , . . . , r } is such that rank ( G H ( S )) = rank ( S ) (12) then ker (cid:0) G H ( S ) T (cid:1) = ker (cid:0) S T (cid:1) (13) and thus in particular the stoichiometric matrices G H ( S ) and S define thesame SCCs. itle Suppressed Due to Excessive Length 9 Proof
From Definition 7 it follows that ker (cid:0) S T (cid:1) ⊆ ker (cid:0) G H ( S ) T (cid:1) . Indeed, if γ ∈ ker (cid:0) S T (cid:1) then = γ T S = (cid:0) γ T S , . . . , γ T S r (cid:1) where S j , j = 1 , . . . , r , denotes the j − th column of S . In particular γ T S j = 0for all j = 1 , . . . , r , j (cid:54)∈ H , that is γ ∈ ker (cid:0) G H ( S ) T (cid:1) .Additionally, equation (12) implies that ker (cid:0) G H ( S ) T (cid:1) and ker (cid:0) S T (cid:1) havethe same dimension and thus are equal. Corollary 1
Consider a CRN described by the system of ODEs (1) and sat-isfying the global stability condition of Definition 5. Let x and y such that SC ( x ) = SC ( y ) . If H ⊆ { , . . . , r } is such that equation (12) holds and theCRN having stoichiometric matrix G H ( S ) satisfies the stability condition thenthe mutated trajectories starting from x and y lead to the same steady state.Remark 7 This holds in particular if x and y are replaced by the initial state x and the corresponding steady state x e , respectively. When considering theCRN identified by G H ( S ) the trajectories starting from x and x e lead to thesame mutated steady state.3.3 Concatenation of mutationMany cancers arise by effect of a series of mutations accumulated in the cellover time. Here we generalize the results discussed in the previous sections tomodel the simultaneous action of multiple mutations on a given cell.Consider for example a cell affected by two mutations resulting in the LoFof the elemental species A j and A j . The most natural approach to quantifythe combined effect of the two mutations is to start from a concentrationvector x e modeling the (steady) state of a cell in physiological condition, andthen apply the procedure described in Section 3.1 for each single mutation oneafter the other. Specifically, first we project x e through P L j and we computethe steady state of the trajectory started from P L j ( x e ). Then we use P L j toproject the novel, mutated steady state, and we compute the trajectory startedfrom such a projection. The obtained trajectory will belong to the SCC (cid:8) y ∈ R n + s.t. Ny = ˜ c (cid:9) (14)where ˜ c has been obtained by setting to zero the j − th and the j − th elementof c := Nx e .The same result could have been obtained by reversing the order of themutations. First we project x e through P j , then we apply P j to the corre-sponding steady state. The resulting trajectory will still belong to the SCCdescribed by equation (14) and thus will lead to the same steady state ob-tained in the previous case, provided that the CRN satisfy the global stabilitycondition.This result easily generalizes to an arbitrary set of combined GoF and LoFmutations through the following definition. Definition 8
Consider an elemented CRN described by the system (1) withstoichiometric matrix S . Consider a set of (cid:96) mutations, resulting in the LoFof the elemental species A j , . . . , A j (cid:96) , and a set of q GoF mutations, resultingin the suppression of the sets of reactions H , . . . , H q ⊆ { , . . . , r } .The combined effect of the considered mutations is quantified by computingthe asymptotically stable state of the system of ODEs (1) with stoichiometricmatrix G H q ◦ · · · ◦ G H ( S ) (15)and initial condition x (0) = P L j(cid:96) ◦ · · · ◦ P L j ( x e ) (16)where ◦ denotes function composition and x e are the values of species concen-tration reached by the cell in physiological condition. Theorem 3
The steady state computed through the procedure described inDefinition 8 does not depend on the order of the composition in the equations(15) and (16) provided that rank (cid:0) G H q ◦ · · · ◦ G H ( S ) (cid:1) = rank( S ) (17) and that all the involved CRNs satisfy the global stability condition of Defini-tion 5.Proof The thesis follows by observing that the equations (15) and (16) do notdepend on the order of the composition. Indeed the mutated stoichiometricmatrix G H q ◦ · · · ◦ G H ( S ) is defined by setting to zero the columns of S corre-sponding to the reactions in H ∪ · · · ∪ H q that clearly does not depend on theorder in which the GoF mutations are considered. Additionally if condition(17) holds then G H q ◦ · · · ◦ G H ( S ) and S define the same SCCs as shown inTheorem 2. The initial condition (16) defines the unique SCC SC ( P L j(cid:96) ◦ · · · ◦ P L j ( x e )) = (cid:8) x ∈ R n + s.t. Nx = ˜ c (cid:9) where ˜ c has been obtained from c := Nx e by setting to zero the j i − th element,for all i = 1 , . . . (cid:96) . A different order of the LoF mutations results in a differentinitial concentration vector on the same SCC. Therefore, according to Theorem1 it leads to the same mutated asymptotically steady state. Remark 8
Very recent experimental studies on cancer cells have shown theimportance of the order of mutations on the development of cancer diseases(Levine et al., 2019). This seems to be in contrast with the results statedin Theorem 3. However, that theorem follows from a model that does notaccount for the selection process induced on the mutated cell by the externalenvironment. Instead, our model refers to the specific G1-S transition of acancerous cell and it mimics the effects of accumulated mutations just for thatphase. itle Suppressed Due to Excessive Length 11
Remark 9
Coherently to the previous Remark, the process described in Defi-nition 8 is equivalent to apply each considered mutation individually one afterthe other. If the hypotheses of Theorem 3 hold, changing the order of themutation affects the covered trajectory but will lead to the same final steadystate. β , WNT and EGF growthfactors. Following activation, a cascade of chemical reactions proceeds. Therelevant output of the system is given by the concentrations of a number oftranscription factors (activators or repressors). CRC-CRN involves 8 constant,and 411 varying chemical species. The constant species model non-consumablechemicals: three constants describe the growth factors. Internal interactionsare modeled by 339 reversible and 172 irreversible reactions, for a total of r = 850 reactions. Chemical kinetics is based on mass action law, so that 850rate constants are considered. The state of CRC-CRN is described by n = 419variables resulting from the concentrations of the 411 internal species, plus8 additional variables, of null time derivative, accounting for the constantinflows. The state variables satisfy a system of 419 polynomial ODEs as inequation (1). In particular the monomials defining the reaction fluxes v ( x , k )are quadratic, since the reactants considered depend at most on two chemicalspecies. The rate constants enter the system as real and positive parameters.4.2 Conservation laws and elemental variables of the CRC-CRNWe have obtained a basis of ker( S T ) consisting of p = 81 semi-positive con-servation vectors, providing an independent set of 81 CLs, all of them being regarded as moiety CLs. In particular, we point out that CRC-CRN satisfiesthe condition p = n − rank( S ). As expected, 8 CLs correspond to constancyof original non-consumable chemical species, so that they can be consideredas trivial. The remaining 73 CLs describe effective properties of the system.It is seen by inspection that the CRC-CRN is weakly elemented, with 9chemical species not involved in the conservation laws. Up to the 8 constantspecies, the list of elemental species coincides with that of the basic speciesof (Tortolina et al., 2015), which were defined as consisting of proteins in freeform, which may bind to other species in order to form derived compounds orsecondary species.A simple, rather typical, example of conservation law within our CRC-CRNisNLK + NLKP + NLKP − TCFLEF + NLKP − Pase10 + TAKP − TAB − NLK = c (18)where c is a given constant; here NLK denotes the concentration of the ele-mental species “Nemo like kinase” (Tortolina et al., 2015); NLPP is the phos-phorylated form of NLK, while an expression like NLKP − TCFLEF refers tothe compound formed by NLKP and TCFLEF. Equation (18) expresses con-servation in time of the total number of molecules of NLK belonging to thecompounds entering the combination in the left side, as discussed in § c ,which is determined by the initial conditions, may be regarded as a counter ofthe conserved molecules of NLK.4.3 Global stability of CRC-CRNThere are a number of results available for equilibrium points and their stabil-ity properties (Chellaboina et al., 2009; Domijan and Kirkilionis, 2008; Fein-berg, 1995; Yu and Craciun, 2018), but they cannot be applied straightfor-wardly to this CRC-CRN, essentially for two reasons. In fact, some of theseresults require very technical hypotheses that cannot be verified in the case ofour ODEs system, due to its dimension and complexity. Further, other resultscannot be applied in the case of open systems, like the one considered in thispaper. Therefore, we will make use of numerical simulations to support thefollowing conjecture. Conjecture 1
CRC-CRN satisfies the global stability condition introduced inDefinition 5.Numerical verification.
We defined 5 different SCCs by randomly selecting thevalues of the total concentrations c . Specifically, each element c j , j = 1 , . . . , p ,has been drawn from a log -uniform distribution on [10 − , ]. For each ofthe obtained SCC we generated 30 initial points x ( k )0 ∈ (cid:8) y ∈ R n + s.t. Ny = c (cid:9) , k = 1 , . . . , itle Suppressed Due to Excessive Length 13 as follows. – First, we dealt with the species that do not belong to any conservation law.The value of their initial concentration were log − uniformly drawn fromthe interval [10 − , ]. – Then we considered all the other chemical species but the elemental species.After randomly permuting their order, for each species i :1. we randomly selected an initial concentration value below the upper-bound imposed by the total concentrations available in the conservationlaws involving it. More in details we set x ( k )0 ,i = u · min j ∈ Γ ( i ) c j γ ji , where u was uniformly drawn from [0 ,
1] and Γ ( i ) = { j ∈ { , . . . , p } s.t. γ ji (cid:54) = 0 } . is the set of all the conservation laws involving the i − th species.2. we updated the total concentration available in each conservation lawby removing the amount already filled by the i − th species, i.e. c j ← c j − γ ji x ( k )0 ,i for all j ∈ Γ ( i ). – Finally we considered the elemental species. By exploiting the fact thateach elemental species belongs to only one conservation law, j , we set x ( k )0 ,i = c j where c j is the value of the total concentration still available after theprevious step.For each of the 30 points generated with the described procedure, we usedthe matlab tool ode15s (Shampine and Reichelt, 1997) to integrate the systemof ODEs (1) on the interval [0 , . · ], with initial condition x (0) = x ( k )0 . Thevalue of the solution at the last time–point was considered as the correspondingasymptotic steady state x ( k ) e .Since the initial values x ( k )0 all belong to the same SCC they should leadto the same steady state. This was verified by computing for each speciesthe coefficient of variation of the equilibrium values across the 30 trajectories.Namely, for each species we computed ε e ( i ) = (cid:113) (cid:80) k =1 ( x ( k ) e,i − µ e ( i )) µ e ( i ) (19)where µ e ( i ) = (cid:80) k =1 x ( k ) e,i .As a comparison, for each species we also computed the coefficient of vari-ation ε ( i ) across the initial values x ( k )0 ,i , k = 1 , . . . , Fig. 1
Average and standard error of the mean over the non–constant species of the coeffi-cient of variation across 30 different initial conditions (upper panel) and the correspondingasymptotically stable states (lower panel). Each bar corresponds to the results obtainedwith a different SCC. Note the different scale on the y axes. Figure 1 shows the averaged coefficient of variations obtained with the 5considered SCCs. Fixed a SCC, the coefficient of variation across the initialvalues ε is always around 2.5, while the coefficient of variation across thesteady states ε e is at least one order of magnitude lower. Moreover, the lattershows higher differences across the SCCs. This is mainly due to the fact thatdepending on the SCC, few species may require a longer time to reach theasymptotically stable state and thus may show a higher variation when thesystem of ODEs is integrated in the fixed time-interval [0 , . · ].4.4 LoF mutation of TBRIIWe investigated the impact on CRC-CNR of a mutation resulting in the LoFof the elemental species TBRII.First, we computed the value of the concentration vector modeling thephysiological steady state of the cell prior to mitosis. To this end, we integratedthe system of ODEs (1) on the interval [0 , . · ], with initial condition x (0) = x , where x is the ideal state defined by setting the initial concentrations ofthe elemental species as in Tortolina et al. (2015). The value of the steadystate in the physiological cell was defined as the value of the solution at thelast time–point. itle Suppressed Due to Excessive Length 15 Fig. 2
Value of the relative difference δ i between the steady states in the physiological celland in the cell affected by LoF mutation of TBRII. Black and red lines are obtained whenthe system of ODEs for the mutated cell is solved with initial condition x (0) = P L j ( x )and x (0) = P L j ( x e ), respectively. By using Definition 6, we then defined the operator P L j associated to theLoF of TBRII. The steady state of the mutated cell was computed by inte-grating the system of ODEs (1) with two different initial conditions, namely x (0) = P L j ( x ) and x (0) = P L j ( x e ). As in the previous step, in both cases wesolved the system on the interval [0 , . · ] and we defined as steady statethe value of the solution computed at the last time–point.According to Theorem 1 the two trajectories, starting from P L j ( x ) and P L j ( x e ), should lead to the same steady state. This result was numericallyverified by computing for both the trajectories the relative difference δ i = x me,i − x e,i x e,i , i = 1 , . . . , n (20)where x me is the steady state in the mutated cell.As shown in Figure 2, the results obtained with the two different initialconditions coincide and thus Theorem 1 is verified. Moreover, Figure 2 showsthe effect of the LoF of TBRII on the concentrations of all the chemical species.Indeed, δ i < i − th speciesin the mutated cell is lower than in the healthy cell. On the contrary δ i > i − th species is higher in the mutated cell.On the other hand, different values of the initial conditions lead to differenttrajectories. In Figure 3 we compared the time required by the two trajectoriesto reach a stable state. To this end for each time point t we computed || ˙ x ( t ) || ∞ = || Sv ( x ( t ) , k ) || ∞ . (21)Since x e is already a steady state for the CRC-CRN, when computing theprojected value P L j ( x e ) the species not affected by the LoF of TBRII maintaintheir stable values. As a consequence, the trajectory starting from P L j ( x e )requires a lower number of iterations to reach a value of || ˙ x ( t ) || ∞ closer to 0. Fig. 3
Infinity norm of the derivative ˙ x of the concentration vector as function of time. Blackand red lines show the results obtained by solving the ODEs system (1) with initial condition x (0) = P L j ( x ) and x (0) = P L j ( x e ), respectively, P L j being the operator associate to theLoF mutation of TBRII G H , where the H is defined to include all the reactions involved in the deactivation of BRAF ∗ .We recall that BRAF ∗ is the activated form of BRAF, consisting in the phos-phorylation of a specific amino acid. Thus it is assumed that the mutated formof BRAF is still subject to activation in BRAF ∗ , while inactivation of BRAF ∗ is blocked.As described in Tortolina et al. (2015), Supplementary materials, the deac-tivation of BRAF ∗ is regulated by the phosphatase Pase1 through the followingset of reactions: BRAF ∗ + Pase1 k f (cid:29) k r BRAF ∗− Pase1BRAF ∗− Pase1 k → BRAF + Pase1where a simplified notation has been adopted for the rate constants. To blocksuch a deactivation process while respecting the condition described by equa-tion (12) in Theorem 2 we removed the two forward reactionsBRAF ∗ + Pase1 k f → BRAF ∗− Pase1 k → BRAF + Pase1 , that is we defined a novel stoichiometric matrix G H ( S ) were the correspondingcolumn of S have been set equal to 0.Similarly to what we have done in the previous section, we computed theasymptotically stable states of the trajectories obtained solving the system ofODEs defined by the stoichiometric matrix G H ( S ) with two different initialcondition, namely x (0) = x and x (0) = x e , x and x e being the ideal andsteady state for the physiological cell.Figure 4 shows that the two trajectories lead to the same steady state asstated in Corollary 1. Additionally, as done in the previous section, in Figure itle Suppressed Due to Excessive Length 17 Fig. 4
Value of the relative difference δ i between the steady states in the physiological celland in the cell affected by GoF mutation of BRAF. Black and red lines are obtained whenthe mutated system of ODEs defined by the stoichiometric matrix G H ( S ) is solved withinitial condition x (0) = x and x (0) = x e , respectively. Fig. 5
Infinity norm of the derivative ˙ x of the concentration vector as function of time whenthe stoichiometric associated to GoF of BRAF is employed in the ODEs system (1). Blackand red lines show the results obtained by setting x (0) = x and x (0) = x e , respectively. S , and thus modifying thesystem of ODEs associated to the CRC-CRN. Although x e was a stable statefor the system of ODEs modeling the cell in physiological condition, for mostof the species the corresponding concentration x e,i is no longer an equilibriumvalue for the new system. For this reason the two trajectories, starting from x e and x , show similar behaviors of || ˙ x ( t ) || ∞ as function of the time t in whichthe solution of the system is computed. In this paper we have considered a system of ODEs modeling a CRN, withemphasis on the mutual relationships between moiety CLs and stoichiometricsurfaces. CLs have led to the characterization of weakly elemented CRN andthe definition of ideal states, which are crucial in the treatment of mutations.Confining attention to systems satisfying the global stability condition, wehave considered two classes of mutations, LoF and GoF, often found in cancer cells. Mutations have been modeled as projection operators and basic prop-erties of mutated networks have been analyzed, aiming at the determinationof mutated equilibrium states as asymptotic steady limits of trajectories. Thenew model for mutations allows for a simple treatment of their combinations,showing that the resulting equilibrium is independent of the order.As an application of the previous results we have investigated a systemof ODEs describing the G1-S phase of a CRC cell. Due to the huge numberof variables involved, the analysis has been based on numerical simulations.We have found 81 independent, linear, moiety CLs; they have been applied tosupport the conjecture that CRC-CRN satisfies the global stability condition.It has also been found that the elemental chemical species coincide with thebasic species defined by Tortolina et al. (2015) on biochemical grounds. Also,a mutation by LoF, and another mutation by GoF have been examined indetail.We are aware that this computational approach to CRN is devised to cap-ture few, although decisive, aspects of the G1-S transition point of a cell,as well as the modifications of the network induced by cancer mutations. Wethink that the model can provide a guide to future experimental research basedon the interpretation of results of simulations, particularly in the case of thedevelopment and optimization of targeted therapies agains already mutatedcells. For example, if the simulation predicts a significant increase of the mu-tated concentration of a specific metabolite with respect to its physiologicalvalue, then that metabolite can be regarded as a potential drug target. Morein general, an analysis of the mutated profile of the simulated CRN in theG1-S phase should provide hints about convenient choices of targets for avail-able drugs, together with a framework for the simulation of their consequences(Facchetti et al., 2012; Torres and Altafini, 2016).Possible development of our model may involve different directions. For ex-ample, we could extend it to examine alterations of mRNA induced by changesfrom physiologic to mutated equilibrium (Castagnino et al., 2016). Further, inthe model we have not considered possible dependence of the literature dataon temperature, pH, or other conditions. Therefore, the scheme is not capa-ble of describing the individual response of a fixed cell; rather, it attempts atsimulating the behavior of a homogeneous set of cells; from this viewpoint,the solution of the system of ODEs may be interpreted as providing an aver-age dynamic of the intra-cellular answer. Also, the scheme should be enrichedin order to account for the impact of the extra-cellular micro-environment,which implies to modify the model in order to account for growth-inducedsolid stresses and pressure, nutrients and oxygen supply, and blood perfusion(Caviglia et al., 2014; Jones and Chapman, 2012; Markert and Vazquez, 2015).As a final comment, we remark that the rate constants in the CRC-CRNconsidered in this paper are assumed as known and estimated from the scien-tific literature (Tortolina et al., 2015). A more reliable determination of theseconstants can be obtained by solving a non-linear ill-posed problem (Berteroand Piana, 2006; Engl et al., 2009) in which the input data are provided by ad hoc conceived biochemical experiments. The setup of such experiments and itle Suppressed Due to Excessive Length 19 the realization of an optimization method for the numerical solution of thisinverse problem will be part of future research.
Acknowledgments
MP has been partially supported by Gruppo Nazionale per il Calcolo Scien-tifico. The authors kindly acknowledge Prof. Silvio Parodi, Prof. Silvia Raveraand Dr. Mara Scussolini for their useful comments and feedback.
References
Anderson MW, Moss JJ, Szalai R, Lane JD (2019) Mathematical modelinghighlights the complex role of AKT in TRAIL-induced apoptosis of colorec-tal carcinoma cells. iScience 12:182–193Bertero M, Piana M (2006) Inverse problems in biomedical imaging: modelingand methods of solution. In: Complex systems in biomedicine, Springer, pp1–33Castagnino N, Maffei M, Tortolina L, Zoppoli G, Piras D, Nencioni A, MoranE, Ballestrero A, Patrone F, Parodi S (2016) Systems medicine in colorec-tal cancer: from a mathematical model toward a new type of clinical trial.WIREs Syst Biol Med 8(4):314–336Caviglia G, Morro A, Pinamonti N (2014) The Klein–Gordon equation inmixture models of tumour growth. Phys Lett A 378(48):3607–3613Chellaboina V, Bhat SP, Haddad WM, Bernstein DS (2009) Modeling andanalysis of mass-action kinetics. IEEE Control Systems Mag 29(4):60–78De Martino A, De Martino D, Mulet R, Pagnani A (2014) Identifying allmoiety conservation laws in genome-scale metabolic networks. PloS one 9(7)Domijan M, Kirkilionis M (2008) Graph theory and qualitative analysis ofreaction networks. Net Heter Media 3(2):295Eduati F, Dold`an-Martelli V, Klinger B, Cokelaer T, Sieber A, Kogera F,Dorel M, Garnett MJ, Bl¨uthgen N, Saez-Rodriguez J (2017) Drug resistancemechanisms in colorectal cancer dissected with cell type–specific dynamiclogic models. Cancer res 77(12):3364–3375Engl HW, Flamm C, K¨ugler P, Lu J, M¨uller S, Schuster P (2009) Inverseproblems in systems biology. Inverse Probl 25(12):123014Facchetti G, Zampieri M, Altafini C (2012) Predicting and characterizing se-lective multiple drug treatments for metabolic diseases and cancer. BMCSyst Biol 6(1):115Feinberg M (1987) Chemical reaction network structure and the stability ofcomplex isothermal reactorsI. The deficiency zero and deficiency one theo-rems. Chem Engin Sci 42(10):2229–2268Feinberg M (1995) The existence and uniqueness of steady states for a classof chemical reaction networks. Arch Rat Mech Anal 132(4):311–370Griffiths AJ, Wessler SR, Lewontin RC, Gelbart WM, Suzuki DT, Miller JH,et al. (2005) An introduction to genetic analysis. Macmillan
Haraldsd´ottir HS, Fleming RM (2016) Identification of conserved moieties inmetabolic networks by graph theoretical analysis of atom transition net-works. PLoS Comput Biol 12(11):e1004999Hochman G, Halevi-Tobias K, Kogan Y, Agur Z (2017) Extracellular inhibitorscan attenuate tumorigenic wnt pathway activity in adenomatous polypo-sis coli mutants: Predictions of a validated mathematical model. PloS one12(7):e0179888Jones GW, Chapman SJ (2012) Modeling growth in biological materials. SiamRev 54(1):52–118Jordan JD, Landau EM, Iyengar R (2000) Signaling networks: the origins ofcellular multitasking. Cell 103(2):193–200Karin M, Smeal T (1992) Control of transcription factors by signal transduc-tion pathways: the beginning of the end. Trends Biochem Sci 17(10):418–422Kohn KW, Aladjem MI, Weinstein JN, Pommier Y (2006) Molecular interac-tion maps of bioregulatory networks: a general rubric for systems biology.Mol Biol Cell 17(1):1–13Kschischo M (2010) A gentle introduction to the thermodynamics of bio-chemical stoichiometric networks in steady state. Eur Phys J Special Topics187(1):255–274Lemieux E, Cagnol S, Beaudry K, Carrier J, Rivard N (2015) Oncogenic KRASsignalling promotes the Wnt/ β -catenin pathway through LRP6 in colorectalcancer. Oncogene 34(38):4914–4927Levine AJ (2019) Targeting therapies for the p53 protein in cancer treatments.Ann Rev Cancer Biol 3:21–34Levine AJ, Jenkins NA, Copeland NG (2019) The roles of initiating truncalmutations in human cancers: the order of mutations and tumor cell typematters. Cancer cell 35(1):10–15Li Y, Zhang Y, Li X, Yi S, Xu J (2019) Gain-of-function mutations: an emerg-ing advantage for cancer biology. Trends Biochem SciMarkert EK, Vazquez A (2015) Mathematical models of cancer metabolism.Cancer & metabolism 3(1):14Roy M, Finley SD (2017) Computational model predicts the effects of targetingcellular metabolism in pancreatic cancer. Front Physiol 8:217Schilling CH, Letscher D, Palsson BØ (2000) Theory for the systemic definitionof metabolic pathways and their use in interpreting metabolic function froma pathway-oriented perspective. J Theor Biol 203(3):229–248Schuster S, H¨ofer T (1991) Determining all extreme semi-positive conservationrelations in chemical reaction systems: a test criterion for conservativity. JChem Soc Faraday Trans 87(16):2561–2566Sever R, Brugge JS (2015) Signal transduction in cancer. Cold Spring HarbPerspect Med 5(4):a006098Shampine LF, Reichelt MW (1997) The MATLAB ODE suite. SIAM J SciComput 18(1):1–22Shinar G, Alon U, Feinberg M (2009) Sensitivity and robustness in chemicalreaction networks. Siam J Appl Math 69(4):977–998 itle Suppressed Due to Excessive Length 21itle Suppressed Due to Excessive Length 21