In-silico modelling of the mitogen-activated protein kinase (MAPK) pathway in colorectal cancer: mutations and targeted therapy

Introduction: Chemical reaction networks (CRNs) are powerful tools for describing the complex nature of cancer’s onset, progression, and therapy. The main reason for their effectiveness is in the fact that these networks can be rather naturally encoded as a dynamical system whose asymptotic solution mimics the proteins' concentration profile at equilibrium.

Methods and Results: This paper relies on a complex CRN previously designed for modeling colorectal cells in their G1-S transition phase and presents a mathematical method to investigate global and local effects triggered on the network by partial and complete mutations occurring mainly in its mitogen-activated protein kinase (MAPK) pathway. Further, this same approach allowed the in-silico modeling and dosage of a multi-target therapeutic intervention that utilizes MAPK as its molecular target.

Discussion: Overall the results shown in this paper demonstrate how the proposed approach can be exploited as a tool for the in-silico comparison and evaluation of different targeted therapies. Future effort will be devoted to refine the model so to incorporate more biologically sound partial mutations and drug combinations.

1 Introduction

For many years, the primary sources of anti-cancer intervention just consisted in chemiotherapy, where specific drugs are used to kill cancer cells, and surgery. However, surgery is not always feasible, while the broad-spectrum of chemioterapy drugs, which attack indiscriminately cancer and fast-growing healthy cells, often results in high toxicity (Lowenthal and Eaton, 1996). To overcome this low specificity of chemiotherapy, novel biology-based approaches have been introduced in routine cancer therapies. In particular, following the advances of human genome sequencing, in the last 20 years an increasing number of molecular targeted drugs have been introduced. This kind of therapeutic intervention aims at slowing down cancer progression and metastasis by targeting specific molecules somehow involved in the genetic alterations that underlie cancer onset (Lee et al., 2018; Bedard et al., 2020; Zhong et al., 2021). Additionally, synergies of multiple targeted drugs combined in a single therapy have been investigated to reduce resistance to single-agent therapies (Jin et al., 2023). Identifying novel molecular targets and optimizing dosage and combination of the corresponding drugs remains a challenging problem, where the number of possible therapies to be tested vastly exceeds clinical resources, in terms of both financial and time resources. In this scenario, systems biology models could play a crucial role in identifying the most promising candidates for clinical trials and in elucidating the molecular mechanisms underlying targeted drugs synergies (Chen et al., 2015; Rocca and Kholodenko, 2021).

At a worldwide level, colorectal cancer (CRC) is the third most frequent cancer in male population and the second one among women (Biller and Schrag, 2021; Sung et al., 2021; Xi and Xu, 2021). Screening programs have contributed in reducing the incidence of later-onset cases while an alarming increase in early-onset cases and in corresponding CRC-related mortality among younger people have been observed (Saad El Din et al., 2020; Sinicrope, 2022; Wu and Lui, 2022). At a molecular level CRCs are highly heterogeneous pathologies with differences across age groups. It has been estimated that five to ten tumor-specific driver mutations usually concur in individual cancers and that the most frequent alterations in CRC pertain to TP53, APC, KRAS, PTEN, SMAD4, PIK3CA, BRAF, and AKT (Tortolina et al., 2015; Tariq and Ghias, 2016; Anderson et al., 2019). Among these genes, KRAS, APC, SMAD4, and TP53 belong to four different pathways, namely MAPK, WNT, TGFβ, and TP53, each one acting at a different functional stage of cell development, ranging from stem cell renewal to cell growth, division and apoptosis (Armaghany et al., 2012). The first targeted therapies for CRCs, namely cetuximab (Jonker et al., 2007), which inhibits the epidermal growth factor receptor (EGFR), and bevacizumab (Los et al., 2007; Rosen et al., 2017), against the vascular endothelial growth factor A (VEGF-A), were approved by the Food and Drug Administration (FDA) in 2004. Since then, various pathways have been proved to offer ideal sites for targeted therapies, and an increasing number of novel agents have been developed (for a comprehensive review we refer to Tiwari et al. (2018); Xie et al. (2020), and references therein). However, to date only a few CRC-related pathways have been successfully inferred due to the complex signaling network that makes it hard to completely inhibit specific biological interactions. As a consequence, many proposed therapies have not passed the preclinical status or the phase I trial, highlighting the need for systems biology models capable of guiding the choice of which drugs to test so as to avoid waste of resources.

It is clear that a mathematical model aiming at capturing the complex nature of CRC onset, progression, and therapy cannot consider altered and targeted proteins and corresponding pathways in isolation but must integrate them within proper chemical reaction networks (CRNs) (Tortolina et al., 2015). An extensive CRN for CRC has been recently introduced for modelling signal transduction during the G1-S transition phase in colorectal cells (Tortolina et al., 2015; Castagnino et al., 2016). Such a CRN, henceforth denoted as CR-CRN, comprises 10 different pathways, including all four previously mentioned ones, for a total of 419 proteins interacting in 850 chemical reactions. Following standard mathematical procedures based on the law of mass action (Feinberg, 1987; Chellaboina et al., 2009; Yu and Craciun, 2018), the CR-CRN has been mapped into a system of 419 autonomous ordinary differential equations (ODEs) whose solutions describe the behaviour of protein concentrations, which evolve in time until the network reaches an equilibrium (Ingalls, 2013; Sommariva et al., 2021a). It has been conjectured that the CR-CRN satisfies the so called global stability condition (Sommariva et al., 2021a), meaning that a globally asymptotically stable state exists (that is also an equilibrium of the system) once fixed the initial values of the protein concentrations or, more precisely, once fixed the total moiety within the conservation laws of the dynamical system (Shinar et al., 2009; De Martino et al., 2014). By exploiting these properties, Sommariva et al. (2021a) proposed a formal mathematical model to incorporate in the system various, possibly concurrent genetic alterations resulting in a loss or gain of function (L/GoF) of some of the proteins in the network. From a mathematical viewpoint both LoF and GoF mutations are modelled as projection operators acting on the initial concentrations and the stoichiometric matrix of the system, respectively, while the function composition of these operators mimics the action of multiple concurrent mutations.

A rich plethora of information can be derived from the analysis of the solutions of the dynamical systems associated to the original CR-CRN and to its mutated forms. For example, feedback effects naturally emerge by comparing the time-courses of the individual protein concentrations or by studying the corresponding time-varying reaction fluxes (Sommariva et al., 2021b). More importantly, local and global effects induced on the network by LoF and GoF mutations can be quantified by computing the relative difference between the protein concentrations at the equilibrium point of the original CR-CRN and that at the equilibrium of the network obtained by applying the proper projection operator(s) (Sommariva et al., 2021a). Additionally, analysis of the sensitivity of protein concentrations at equilibrium with respect to the values of the kinetic parameters of the dynamical model may help in identifying the specific reaction or subnetwork mostly affected by each genetic alteration (Biddau et al., 2023).

The proposed model has been used to simulate the functional alterations induced on the CR-CRN by some of the mutations more commonly found in CRC, including the GoF of genes PI3K, KRAS, and BRAF, and the LoF of genes PTEN, AKT, and TGFβRII (Sommariva et al., 2021a; Sommariva et al., 2021b; Biddau et al., 2023), and a first attempt in modelling the action of Dabrafenib, a drug targeting BRAF, has been performed by Sommariva et al. (2021b). The obtained results have been extensively compared with results previously published in literature. However, in all those studies only mutations resulting in a complete LoF or in the highest possible value of GoF of the corresponding protein have been considered. In the present paper we applied the proposed CR-CRN in two novel scenarios: (i) the simulation of global and local effects of different genetic mutations resulting in different levels of alteration of the functional activity of the corresponding proteins; (ii) the in silico modelling and dosage of a multi-agent targeted therapy. Toward this end, we focused on mutations and drugs involving proteins directly or indirectly related to the mitogen-activated protein kinase (MAPK), whose overexpression plays an important role in CRC-progression (Fang and Richardson, 2005; Horst et al., 2012; Baharudin et al., 2017; Stefani et al., 2021). More specifically, we quantified the global effects induced on the whole CR-CRN and the local effects induced on the molecules of MAPK by the LoF of PTEN, various levels of GoF of KRAS, and their combination. Finally we investigated the synergies between Dabrafenib and Trametinib, a combination therapy that has demonstrated good results both in terms of progression free survival and response rate (Xie et al., 2020).

2 Materials and methods

The equilibrium concentration of the chemical species in both physiological and diseased conditions is interpreted as an asymptotically stable state of the dynamical system

$\left\{\begin{array}{c}\dot{\mathbf{x}}\left(t\right)=\mathbf{S}\mathbf{v}\left(\mathbf{x}\left(t\right),\mathbf{k}\right)\hfill \\ \mathbf{x}\left(0\right)={\mathbf{x}}_0,\hfill \end{array}\right.(1)$

which is obtained by applying the mass action kinetics to the network represented in the molecular interaction map (see, for example, Kohn, 1999; Pommier et al., 2004; Krogan et al., 2015; Kondratova et al., 2018; Broyde et al., 2021). In Equation 1, $\mathbf{x}(t)={(x_1(t),\dots ,x_n(t))}^T$ is the vector whose components are the concentrations of the n proteins contained in the network; $\mathbf{k}={(k_1,\dots ,k_r)}^T$ is the vector whose components are the rate constants of the r chemical reactions; S is the stoichiometric matrix; v is the vector of reaction fluxes; and x₀ is the auxiliary initial condition.

The first objective of an analysis of the cancer signalling network performed by means of the dynamical system (1) is the characterization of the equilibrium states of the network as the asymptotic behaviour of the solution x = x(t). The main conceptual issue in this respect is that system (1) does not necessarily have a unique equilibrium solution (Feinberg, 1987; Conradi and Flockerzi, 2012; Conradi and Mincheva, 2014; Yu and Craciun, 2018). In order to identify formal assumptions that imply this uniqueness the following process should be considered:

1. Given a solution x = x(t) of (1), a set of p semi-positive conservation vectors satisfying p conservation laws can be identified, which belongs to the kernel of the transpose of the stoichiometric matrix S. The transposed forms of the conservation vectors are used to generate the conservation matrix N.

2. The conservation matrix is said weakly elemented if it contains at least one square submatrix equal to the identity matrix of order p. If this holds, the solution vector x and the conservation matrix N can be re-ordered in such a way that the identity matrix acts on the first p elements of the re-ordered solution vector.

3. If N is weakly elemented, for each initial condition x₀ it is possible to construct p hyperplanes in ${\mathbb{R}}^n$, each one defined by the corresponding conservation law. The stoichiometric compatibility class of x₀ is the intersection of the hyperplanes.

4. As said, in general there is no guarantee that the asymptotic solution of (1) is the same for any ${\mathbf{x}}_0\in {\mathbb{R}}^n$. However, it is possible to formulate a conjecture stating that the asymptotic solution is unique for all initial conditions belonging to the same stoichiometric compatibility class. More precisely, a dynamical system and the corresponding CRN are said to satisfy the global stability condition if such conjecture holds true, i.e. if for every stoichiometric compatibility class there exists a unique globally asymptotically stable equilibrium solution x_e.

In order to make a globally stable CRN supportive for the construction of an in silico model of a colorectal cancer cell, three computational issues should be addressed. First, the dynamical system can be modified in order to implement the presence of single DNA mutations and to compute their impact on the resulting proteomic profile. For example, a LoF mutation, which results in the reduction or even the cancellation of the function of a specific protein, is implemented by projecting the initial concentration values describing the physiological cell onto a new initial state in which the concentrations of the mutated protein and the corresponding compounds are set to zero. On the other hand, a GoF mutation, which enhances the expression of a specific protein, is implemented by setting equal to zero the rate constants corresponding to reactions of de-activation for that protein (Sommariva et al., 2021b).

The second issue is concerned with the numerical computation of the (unique) asymptotic solution of (1). Of course, this can be done by applying a numerical method for the solution of the Ordinary Differential Equations (ODEs) contained in the Cauchy’s problem. However, the dynamical computation of this high-dimensional system is numerically demanding and can be effectively replaced by numerical optimization. In fact, it can be shown that the stationary solution of the Cauchy problem can be determined by an algebraic system, which, in turn, is equivalent to a non-linear root-finding problem. Therefore, in this second approach, iterative schemes can be applied to directly compute the equilibrium state, without the need to approximate the solution of the Cauchy’s problem at each time point. In all simulations presented here, the steady states of the considered CRNs have been computed through a recently developed algorithm called NLPC (non linearly projected combined) method (Berra et al., 2022). NLPC is a constrained root-finding method that combines the Newton approach and the gradient descent direction so that convergence results can be analytically proven. Furthermore, in a large set of simulations NLPC has shown both higher accuracy and higher speed than the classical approach based on the integration of the ODEs system.

Finally, it is important to observe that, on the one hand, system (1) is made of a high number of equations, kinetic parameters, and unknown concentrations; but, on the other hand, that just a smaller subset of such equations is representative of chemical reactions that are actually affected by the cancer somatic mutations. Sensitivity analysis is the technical tool that allows the quantitative assessment of the impact of the kinetic parameters’ uncertainty on the CRN equilibrium state. In its local version, sensitivity analysis is an indicator of the impact that a specific mutation has on the expression of the corresponding protein, and, even more importantly, it is able to identify the sub-networks that are mostly affected by the mutation (Biddau et al., 2023).

Notation. Following the current HUGO Gene Nomenclature Committee (HGNC) guidelines, in this paper we denote gene names in italic, while we used non-italicised gene symbols for the corresponding proteins (Bruford et al., 2020). Some exceptions due to different commonly used nomenclatures are proteins K-Ras¹ and B-Raf². Additionally, throughout the paper, expressions such as GoF of KRAS (LoF of PTEN) denote a mutation of gene KRAS (PTEN) resulting in the gain (loss) of function of the protein K-Ras (PTEN) it encodes.

3 Results

The mathematical model described in Sommariva et al. (2021a) has been applied by Sommariva et al. (2021b) to compute modifications in the equilibrium concentrations of the CRC network induced by mutations in a few genes that are rather common in CRC cancerogenesis. In the present paper we focused on a quantitative analysis of the impact of mutations in KRAS, a widely expressed GTP/GDP-binding protein, whose mutated version is found in more than 30% of CRCs. Specifically, in this section we

•Computed the impact of complete and partial mutations of KRAS on the global proteomic profile of the CRC network.

•Compared the modifications induced by the mutated KRAS on the equilibrium with respect to the ones induced by PTEN, a dual protein/lipid phosphatase that triggers the PI3K/PTEN/AKT signalling pathway.

•Investigated to what extent a mutated KRAS impacts on the expression of specific proteins in MAPK.

•Studied both local and global effects of the combination of mutations in KRAS and PTEN.

Then, in section 3.3, the MAPK pathway will be studied as a molecular target for two inhibitors of B-Raf and MEK, respectively.

3.1 Global effects induced by mutations in colorectal cancer

Mutations of KRAS are very common in CRC and belong to the pathway of the main sequence K-Ras/B-Raf/MEK/ERK. In particular, a GoF of KRAS is realized by a modification of the stoichiometric matrix in (1), which, from a chemical viewpoint, corresponds to removing from the CRN all reactions involved in the de-activation of that protein. A complete list of these reactions can be found in Supplementary Appendix S2. An analysis of global effects of the GoF of KRAS on the proteomic profile of a colorectal cancer cell has been studied by Sommariva et al. (2021b) via the computation of

${\delta }_i=\frac{{\tilde{x}}_i^e-x_i^e}{x_i^e},(2)$

where ${\tilde{x}}_i^e$ and $x_i^e$ are the mutated and the physiological equilibrium, respectively. Since δ_i and the difference ${\tilde{x}}_i^e-x_i^e$ have the same sign, the concentration of the i − th protein in the mutated network is either increased if δ_i > 0 or reduced if δ_i < 0. In particular, a value of δ_i equal to −1 means that the function of the i − protein is completely stopped. In more general terms, the value of δ_i quantifies the relative change of the protein concentration, normalized by its value in the physiological network, and thus enables identifying which proteins are more sensitive to the considered mutation. Finally, it is worth noticing that, when defining δ_i as in Equation 2, we exploited that at the physiological equilibrium all the species in the CR-CRN have a non null concentration. Different metrics should be introduced in scenarios where some of the involved complexes tend to zero.

Using this same technique, here we studied the consequence of partial GoF mutations of KRAS, i.e., of a common proportional reduction of the values of the rates of those reactions that model the inhibition of the active forms of K-Ras. Precisely, in the top panel of Figure 1 the three profiles correspond to values of the rates set to 0%, 30%, and 60% of the corresponding physiological values, while the second panel from the top shows how the whole network is affected by a complete mutation of KRAS. In the third panel we compared this latter profile with the one associated to another frequent mutation in CRC, i.e., the LoF of PTEN, which belongs to the distinct PI3K/PTEN/AKT pathway. Finally, in the bottom panel we studied the effect of the combination of the two mutations on the CRN global equilibrium. From a computational point of view, this last simulation has been carried out by modifying the CR-CRN so as to account for both mutations simultaneously and then by computing the steady state of the obtained modified network. Indeed, in a previous work (Sommariva et al., 2021a) we have shown that this approach produces the same equilibrium point that would have been obtained if we had included the two mutations sequentially, regardless of their order. This is because our model does not currently take into account e.g. selection mechanisms induced by the external environment.

FIGURE 1

FIGURE 1. Global effects induced on the CR-CRN by the GoF of KRAS (A,B), the LoF of PTEN (C), and their combination (D). The effects on each protein, i = 1, …, n, of the CR-CRN is quantified by the relative difference δ_i between the concentrations at equilibrium of the physiological and the mutated networks. (A): results for three mutations corresponding to different levels of GoF of KRAS, obtained by setting the rate constants of the reactions inhibiting the active form of K-Ras to 0%, 30%, and 60% of their physiological values. (B,C): global effects of the complete GoF of KRAS and LoF of PTEN, respectively. (D): global effects induced by the combination of the two complete mutations.

The results of Figure 1 show that the whole network is affected by the mutations of KRAS, PTEN, and their combination; that, in particular, significant changes in concentration may involve proteins far from the mutated ones in the graph of the CRN; and that the global equilibrium profile changes smoothly with respect to a smooth variation of the reactions’ rate. Further, on the one hand, the higher impact of the GoF of KRAS, which is visible by comparing the second, third, and fourth panels, is possibly related to the fact that this protein is upstream in the global network. On the other hand, the same three panels seem to show that no intuitive superposition principle applies to the profile associated to the combination of mutations. This is probably a consequence of the fact that the effects of the combination follow from the solution of a combined non-linear problem with appropriate initial conditions. Further, the profile of δ_i in panel (A) implies a growing difference from the physiological equilibrium, as the values of the mutated coefficients tend to zero, that is, as the mutation is intensified.

With reference to global results concerning the mutation of KRAS, Table 1 shows the list of the 10 molecular species exhibiting the highest variations between physiological and mutated equilibrium values. With the only exception of CDC25C, all of them belong to the mitogen-activated protein kinase (MAPK) signalling pathway, consisting of the main sequence K-Ras/B-Raf/MEK/ERK, while Pase3 is a phosphatase acting on phosphorylated ERK proteins. This result, ultimately obtained from the simulation of the mutation of KRAS, agrees with well known aspects of the physiology of the MAPK signaling pathway. Essentially, the activated form of K-Ras is responsible for the transduction of signals, received at the cell surface, to the inside of the cell; this operation is crucial for cell proliferation, growth and differentiation (Morkel et al., 2015; Pappalardo et al., 2016; Porru et al., 2018; Guo et al., 2020; Lavoie et al., 2020). Under the GoF mutation the molecules of K-Ras persist in their active form, which implies an aberrant activation of downstream effectors as B-Raf, MEK, ERK, and leads to a malignant behaviour of the cell. Thus the whole signaling transduction pathway MAPK is uncontrollably triggered by the alteration of KRAS and may lead to out of control cell proliferation. The protein ERK exhibits the maximum difference between mutated and physiological equilibrium values; further, its active (doubly phosphorylated) form p-p-ERK, which occupies the third place in the list of Table 1, is a well known “master regulator of cell behaviour, life and fate” (Lavoie et al., 2020), being deeply involved in cellular responses as cell proliferation, survival, growth, metabolism, migration and differentiation (Guo et al., 2020; Lavoie et al., 2020; Sugiura et al., 2021). This makes the MAPK path a very natural target of drugs which contrast the negative effects of the mutation of KRAS. In particular, the response of active ERK to the delivery of drugs has received special attention (Pappalardo et al., 2016; Santini et al., 2019; Lavoie et al., 2020; Hamis et al., 2021). Finally, we observe that the high value of the relative difference δ_i of the protein p-p-ERK is due to the smallness of the related physiological equilibrium value.

TABLE 1

TABLE 1. Proteins showing the most significant variation in concentration (reported in decreasing order) when the network is affected by the complete GoF of KRAS. Such a change is quantified by the difference between the concentration values at equilibrium in the physiological status and after the mutation (second column), and by the value of the corresponding δ_i (third column).

As a final comment concerning Figure 1 and Table 1, we remark that all these global results follow from the simulated behaviour of the equilibrium solutions of a system of ODEs. Relevant knowledge from cell biology and biochemistry has been applied in the formulation of the basic physiological model of the network and in the description of mutations (Tortolina et al., 2015). Explicit information on the biological consequences induced by mutations and drugs can be recovered by expert users when looking at the local scale, i.e., focusing on the changes induced on the concentration of specific proteins. As an illustrative example, in the next section we present a study conducted on the effect induced by the GoF of KRAS on the proteins of the pathway MAPK.

3.2 Local effects induced by mutations in colorectal cancer

An analysis of local effects of mutations may lead to a more direct connection with the physiology of the signaling network. A few items have already been examined in Sommariva et al. (2021b), mainly focused on the study of the effects of various single-gene mutations on the concentration values of TP53. Here we present a further set of complementary results on the MAPK pathway.

Figure 2 is devoted to a quantitative analysis of the changes of the MAPK cascade induced by mutations. We have dedicated a panel to each of the molecular species K-Ras, B-Raf, MEK, ERK, and to the corresponding active forms K-Ras₋GTP, p-B-Raf, p-p-MEK, p-p-ERK. The histograms inside each panel provide the concentration values at equilibrium in the following conditions: physiological (purple), mutated by GoF of KRAS (blue), mutated by LoF of PTEN (green), combination of the last two (yellow). Note that the range of the y-axis depends on the panel. The consequences of the mutations outlined by the histograms agree with a number of remarks on the physiology of the MAPK cascade that are scattered over the literature (Morkel et al., 2015; Santini et al., 2019; Guo et al., 2020; Lavoie et al., 2020).

FIGURE 2

FIGURE 2. Values of the concentration at equilibrium of the four species belonging to MAPK pathway, namely K-Ras (A), B-Raf (C), MEK (E), ERK (G), and of their active forms, namely K-Ras_GTP (B), p-B-Raf (D), p-p-MEK (F), p-p-ERK (H). The computed equilibria are related to the physiological network (purple) and to the same one but affected by three different mutations: GoF of KRAS (blue), LoF of PTEN (green) and the concurrence of the two (yellow). For ease of visualization a different scale is used on the y-axis of each label.

For example, the first two columns of each panel compare physiological equilibria (purple) with those mutated by the GoF of KRAS (blue). We observe that the mutated value of K-Ras remains rather small, although increased as a consequence of the mutation. As an expected result of the GoF, the value of the active form K-Ras₋GTP is raised to about 70 nM. Also, the physiological concentration of B-Raf is more than halved, while the mutated value of the active form is about five times the physiological one. Similarly, the inactive form of MEK is heavily decreased by the mutation, while the active (doubly phosphorylated) form p-p-MEK is highly augmented; the same remark holds for ERK and p-p-ERK. To summarize, the concentrations of the inactive forms of the key elements of the MAPK pathway are reduced by the mutation of KRAS (with the only exception of K-Ras), while the concentrations of the active forms are heavily enhanced. Consequently, the whole path is abnormally active.

Consideration on the first (purple) and the third (green) columns provides the impact of the LoF of PTEN on the MAPK pathway. It is found that the equilibrium changes of K-Ras, K-Ras₋GTP, p-p-MEK, ERK, p-p-ERK can be overlooked, while small changes can be seen in B-Raf, p-B-Raf, MEK. Since the mutated active forms of the basic elements of MAPK are almost vanishing, we may conclude that there are no sensible effects of the LoF of PTEN on MAPK. An inspection of the fourth columns (yellow) in comparison with the second ones (blue) shows that the combination GoF of KRAS + LoF of PTEN produces essentially the same effects as the mutation of KRAS alone. This is a further indication that the MAPK path is not sensitive to LoF mutations of PTEN.

The marked increase of the concentration of the active form of ERK under the GoF of KRAS, with the consistent reduction of the inactive ERK, leads in particular to abnormal cell proliferation (Lavoie et al., 2020; Hamis et al., 2021), as already observed while commenting Table 1. Further, it is known (Sugiura et al., 2021) that ERK promotes the apoptosis of the cell, and that this property is inhibited by the activation of the species itself. Looking at the histograms, we can observe that the quantity at equilibrium of p-p-ERK grows steeply under the mutation of KRAS, while the concentration level of ERK decreases significantly. According to Yamamoto et al. (2006) and a suggestion by Hamis et al. (2021), this maintained activation of ERK induces downregulation of antiproliferative genes, thus confirming the dangerous character of such a mutation of KRAS.

3.3 Drugs and drug combination targeting the MAPK pathway

The main objective of this section is to simulate the impact of drugs having the MAPK pathway as their target. Our investigation follows the previous analysis of mutations, by considering first the global effects on the whole network, and next the local effects of the MAPK pathway. As to the global effects, we examined the equilibrium states; as to the local effects, we focused on the analysis of the time course of the activated protein p-p-ERK, which is the key factor to assess the main effects arising from the MAPK pathway. In the present approach, the optimal concentration of a drug is determined by the degree of compliance between the equilibrium features of the model modified by the drug, and the corresponding features of the original physiological model. However, the flexibility of the proposed model could be exploited to define other optimal criteria, looking e.g. to the amount of drug that minimizes (even below the level reached in the physiological network) the concentration of proteins coded by certain oncogenes.

To describe the action of a given drug, we enlarge the set of the chemical reactions of the CR-CRN, and the corresponding dynamical system, in order to account for the reactions between the drug and the target molecules. As an immediate consequence, the drug and the associated composites are regarded as additional unknowns. Borrowing from a rather well established literature, we have simulated the action of two drugs: Dabrafenib (DBF) and Trametinib (TMT), with target B-Raf and MEK. As far as reactions are concerned, DBF is modeled as a competitive inhibitor of Raf, while TMT is an allosteric inhibitor of MEK (Morkel et al., 2015; Puszkiel et al., 2019; Sommariva et al., 2021b; Hamis et al., 2021) (more details are given in the Supplementary Appendix S1). The key parameters of the enlarged system are represented by the initial values of the drug concentrations, say c_D for DBF, and c_T for TMT, although also the additional rate constants have to be fixed; the initial values c_D and c_T denote also the total amounts of drug available for the network.

For each inhibitor we examined global and local effects; also, effects of combinations of DBF and TMT are illustrated. In the first set of experiments aimed at investigating global effects induced by these drugs, we assumed they were given concurrently. Indeed, some preliminary tests not shown here suggest that in the current version of our model other approaches of drug administration, e.g., inserting TMT after some time that DBF is in action, does not sensibly change the final equilibrium state reached by the system. In the present approach we considered the mutated equilibrium state as the initial state of the dynamical system modified by the addition of a drug, and we determined the resulting new equilibrium state $x_i^d$. In particular, we found the most appropriate concentrations, c_D and/or c_T, to obtain the closest equlibrium $x_i^d$ closest to the physiological equilibrium $x_i^e$. Figure 3 provides a synthetic description of the changes induced by the action of DBF on the network subject to a GoF mutation of KRAS. This figure has been obtained through the following steps.

•Consider the augmented system formed by the mutated dynamical system, enlarged by addition of the reactions expressing the action of DBF.

•Consider the mutated equilibrium state ${\tilde{x}}_i^e$ and the initial value for the drug concentration c_D as the initial values for the dynamics of the augmented system.

•Determine the equilibrium state $x_i^d$ of the augmented system.

•Compare $x_i^d$ with the physiological equilibrium $x_i^e$ by considering the corresponding relative difference of the concentrations, which here is denoted as

$d_i=\frac{x_i^d-x_i^e}{x_i^e},(3)$

for convenience.

•Investigate the dependence of d_i on c_D, and plot the results.

FIGURE 3

FIGURE 3. Dosage of Dabrafenib (DBF). (A): modified geometric mean G computed on vector $(\mathbf{d})={(d_i)}_i$ of the relative differences between the proteins’ concentration at the equilibrium of the physiological network and in the one obtained by simulating the action of DBF against the GoF mutation of KRAS. Several values for the initial concentration c_D of DBF have been tested and the colored marks corresponds to those further explored in the right panel. (B): Relative difference d_i as function of the proteins quantifying the effect of the GoF mutation of KRAS (black line) and the effect of four initial values of DBF concentration (colored lines).

Clearly, the procedure also applies to the description of changes induced by TMT, with initial value c_T, and to the combination of the two drugs.

In order to determine which is the best quantity of drug to be administered the information in vector d = (d)_i is summarized in a unique numeric value through the indicator

$G\left(\mathbf{d}\right)=\sqrt[n]{\prod _{i=1}^n\left(d_i+10^{-6}\right)}-10^{-6}(4)$

that consists in a modified geometric mean allowing the vector d to contain elements equal to zero (De La Cruz and Kreft, 2018). Such an index is exploited in Figures 3–5 for showing the performance of single drugs DBF and TMT and their combination in dependence of their initial concentration values inside the network. The lower the value of G(d) the closer the equilibrium of the network including the drug(s) is to the physiological equilibrium. The result of the procedure is that the optimal drugs dosage could be determined by minimizing G(d) that is, we assume that the drugs work at their best if the corresponding equilibrium state is close to the physiological (healthy) state. Other metrics could be introduced, e.g., in order to allow the network to surpass the physiological state if there is evidence that this may reduce oncogenic activity.

FIGURE 4

FIGURE 4. Dosage of Trametinib (TMT). (A,B) have been designed as in Figure 3 but considering TMT instead of DBF.

FIGURE 5

FIGURE 5. Dosage of the combination therapy including Dabrafenib (DBF) and Trametinib (TMT). (A): 2-D images representing the modified geometric mean G of the vector $\mathit{d}={(d_i)}_i$ of the relative differences between the proteins’ concentration at the equilibrium of the physiological network and in the network obtained by simulating the concurrent action of DBF and TMT against the GoF mutation of KRAS. Several pairs of values (c_T, c_D) for the initial concentrations of the two drugs have been tested and the colored marks correspond to those further explored in the right panel. (B): Relative difference δ_i as function of the proteins quantifying the effect of the GoF mutation of KRAS (black line) and the effect of four pairs of initial concentration values of TMT, c_T, and DBF, c_D (colored lines).

Figure 3, panel (A) shows the plot of the modified geometric mean G(d) as a function of the initial concentration of DBF c_D, which varies in the interval [0, 100] nM. Therefore, the plot provides an estimate of the drug potential of restoring the healthy state of a cell. The minimum of G corresponds to the initial drug concentration that assures the best healing effect on the mutated network. Figure 3 panel (B) shows a few complete profiles of the relative differences d_i. The black line, which has been inserted for ease of comparison, represents the relative difference between the concentrations at equilibrium under a GoF of KRAS and the physiological equilibrium; notice that the values of δ_i are sorted in decreasing order, to give evidence to the changes induced by the drug. The different colors of the other lines correspond to the profiles of d_i determined by different values of c_D.

The best choice of c_D corresponds to the line closest to the horizontal axis, i.e., c_D = 43.5 nM. This is in essential agreement with the result of panel (A). Notice that the initial values of the concentration of DBF and the values of the rate constants of the related reaction have been changed with respect to those ones used in Sommariva et al. (2021b), as observed in the Supplementary Appendix S1; this explains the slight difference with the optimal value of the initial concentration found in that paper.

In correspondence with the value c_D = 43.5 nM, the concentrations of most proteins involved in network are very close to the values at the physiological equilibrium. We point out a few exceptions: B-Raf, whose concentration is reduced in that its function has been inhibited by the drug; a group of complexes that involve the activated form of K-Ras, that is still overexpressed; the complexes that are products of the reactions removed to simulate the GoF of K-Ras, whose function is thus stopped.

Figure 4 describes the same analysis, this time performed in the case of the effects of TMT on the GoF of KRAS. In this case, the initial concentration values are allowed to vary over a rather large interval, from 200 to 2000 nM, the system being poorly sensitive to variations of c_T. A minimum of G(d) is found at c_T = 1080 nM. Figure 4, panel (B), agrees with the previous result.

Finally, Figure 5 analyzes the effects on the relative error d_i of various combination therapies involving DBF and TMT. The heatmap of panel (A) shows the response to a set of combination therapies. On the horizontal axis, the concentration of TMT varies from 0 to 2000 nM; on the vertical axis, the concentration of DBF describes the interval from 0 to 100 nM. The color of the heatmap reflects the value of G(d) as described by the color map alongside. Specifically, low values of G(d) appear in yellow while high values are in blue. The best result is obtained for c_D ≈ 40 nM and c_T ≈ 240 nM. Thus, a combination of smaller doses of both drugs produces almost the same effects of either a single infusion of DBF, or TMT, delivered at much higher dose.

Panel (B) shows four examples of d_i distribution, corresponding to different choices of the initial conditions c_T and c_D for TMT and DBF. As in the previous analogous representations of the distributions associated with drugs, we have reported by a black line the relative difference δ_i between equilibrium values of the network subject to the GoF of KRAS and the physiological values. Comparison of the results shows that there is a significant convergence between the conclusions drawn from panels (A) and (B) as to the most convenient combination.

The last part of this section is devoted to local considerations on the time course of p-p-ERK. For the ease of reference, we recall that ERK is an elemental conserved variable, according to Sommariva et al. (2021a). Thus, we denote by ERK_tot the corresponding conserved value, and we call activated fraction of ERK the ratio p-p-ERK(t)/ERK_tot. Figure 6 shows the time course of the activated fraction of ERK under the action of DBF in panel (A), TMT in panel (B), and the combination DBF + TMT in panel (C). The initial conditions for the drugs are chosen as c_D ∈ {12.5, 25, 37.5, 50} (nM), c_T ∈ {50, 100, 150, 200} (nM), and the same values are considered for the drug combinations.

FIGURE 6

FIGURE 6. Time-varying behaviour of the concentration of p-p-ERK, active (doubly phosphorylated) form of ERK, in the CR-CRN affected by GoF of KRAS under the action of DBF (A), TMT (B) and their combination (C). In all panels, the activated fraction of ERK is represented. The black dashed line indicates the value of active ERK at equilibrium in physiological conditions, while colored lines refer to four different therapeutic dosage as specified in the legend of each panel where c_D and c_T indicate the initial concentration of DBF and TMT, respectively. A logarithmic scale is used on the time axis.

The graphs of the three panels have a rather similar behaviour. A common feature is that the huge amount of the activated fraction of ERK decreases with time, because of the indirect action of the drugs that target Raf or MEK or both. In the case that only small quantities of delivered drug are available, each panel contains a curve showing an initial slight decrease of the ratio, until it tends to a non vanishing constant value for growing t. On the contrary, the curves associated with higher amounts of drug show an almost horizontal behaviour for a brief initial time interval, after which they decrease very steeply until a rather small value of the activated fraction of ERK is found. Actually, the almost stationary low value is reached in about 13 min or more for DBF and TMT + DBF, and in about 6 min under the action of TMT alone, perhaps because the target MEK of TMT is closer to ERK in the network topology.

4 Discussion

A mathematical model has been recently introduced (Tortolina et al., 2015; Sommariva et al., 2021a) that simulates the behaviour of a CRN describing the information flow inside a CRC cell at the G1-S transition point; furthermore, general procedures to change the model in order to account for GoF and LoF mutations have been proposed and investigated (Sommariva et al., 2021b). In other words, the transition from a healthy to a mutated, cancerous, signaling network has received an appropriate mathematical formulation.

In the present study the mathematical model has been applied in order to analyze the reaction of the network to a GoF mutation of KRAS, a LoF of PTEN, the combination of the two mutations, and a partial GoF mutation of KRAS. The interest toward mutations of KRAS and PTEN comes from the observation that they have been found in approximately 40% and 34% of all CRC cases, respectively (Salvatore et al., 2019; Zhu et al., 2021). A global analysis of pertinent equilibrium states has pointed out the MAPK pathway, with its K-Ras/B-Raf/MEK/ERK cascade, as the target of the main changes of the network, induced by the mutation of KRAS. A corresponding local analysis has provided a suggestive graphical representation of quantitative aspects of the products of activation/inactivation reactions. In particular, the mathematical analysis identifies the activated form of ERK as the focus of changes, in parallel with well known results coming from the biologic side (Guo et al., 2020; Lavoie et al., 2020; Sugiura et al., 2021).

Thus, the species of the MAPK pathway represent the natural target of drugs that tend to contrast the negative consequences of KRAS mutations. Here we have examined the response of the mutated network to administration of DBF and TMT, with respective targets B-Raf and MEK; more precisely, DBF is modeled as a competitive inhibitor of B-Raf, and TMT as an allosteric inhibitor of MEK. As a further related development, the effects of combination of the two drugs at variable doses have also been simulated. In fact, it is well known that drug combinations may counterbalance, e.g., the onset of drug-resistant tumour subclones (Morkel et al., 2015; Santini et al., 2019; Hamis et al., 2021).

The main novelty of the mathematical scheme for the simulation of drug effects is given by the application of three different models, namely, the model for the healthy, the mutated, and the drug loaded network. Here, the mutated network has played a fundamental role. In our approach, we have considered the equilibrium states pertaining to each model. The mutated equilibrium has provided the initial value of the drug loaded model, leading to the associated equilibrium. This last, in turn, has been compared with the healthy equilibrium. A drug, or a drug combination, has been regarded as effective if its drug loaded equilibrium is close to healthy equilibrium, which means that the relative difference d_i is globally small.

Next, a global index has been introduced, dependent on d_i, which selects the amount of drug generating the closest model to physiological equilibrium. The procedure has been applied to assess the optimal concentrations of DBF and TMT. The results have been confirmed by graphical representations of the distribution of the relative differences between drug loaded and physiological equilibrium. Similarly, an optimal drug combination of DBF and TMT has been assessed on the basis of the values of the global index that have been reported on a heatmap. Again, the optimal choice has been validated by the representation of the distributions of relative differences. In details, Figures 3–5 suggest that in the current version of our model there is no combination of the two drugs that allows the system to reach an equilibrium state closer to the physiological one than the one reached by the system when perturbed by either drug alone. However, the combination of the two drugs allows to reduce the optimal dosage, especially the one for TMT, which decreases from about 1080 nM to 240 nM.

At the local level of analysis, we have examined the time course of the ratio p-p-ERK(t)/ERK_tot, representing a fractional measure of the activated ERK. Unlike recent computational approaches as (Pappalardo et al., 2016; Hamis et al., 2021), which, however, consider mutations of KRAS, we have obtained a realistic behaviour of the activated fraction of ERK, whereby the ratio decreases with time, because of drug action. In our opinion this result follows ultimately from the choice of the non-vanishing initial value of p-p-ERK, coincident with that of the mutated equilibrium, and the global effects of the network, incorporating also, e.g., feedback effects.

To further ascertain the realistic behaviour of active fraction of ERK that follows from our approach, we have investigated the changes induced under the assumption that DBF is subject to degradation, while TMT is absent. Precisely, we have considered for DBF a degradation rate equal to 5.79 ⋅ 10^–6 s⁻¹ (Anderson et al., 2019), while all other conditions of the model have been left unaltered. Figure 7 shows the results: the active local fraction of ERK shows an initial value of about 0.5, which is reduced to a very small value in the first hour; the latter is maintained for a rather long time interval, until it reaches again the value 0.5, following the decrease of the drug concentration caused by degradation. If TMT is added to the model, the effects of the degradation of DBF are mitigated. Indeed, as shown by Figure 8, also in this case after about 4800 min the concentration of p-p-ERK started increasing but the ratio $\frac{\text{p}-\text{p}-\text{ERK}}{\mathrm{E}\mathrm{R}{\mathrm{K}}_{\mathrm{t}\mathrm{o}\mathrm{t}}}$ remains below 4 ⋅ 10^–3 nM. Figure 8 also shows that the time at which TMT is administered does not change the final equilibrium state, but impacts on the overall dynamics followed by the system. As an example, if the two drugs are given concurrently, p-p-ERK reached the physiological level faster. Future work will be devoted to a more in-depth study on the capability of our model in capturing the effect of drugs sequencing and/or metromomic therapies.

FIGURE 7

FIGURE 7. Time-varying behaviour of the concentration of p-p-ERK, active (doubly phosphorylated) form of ERK, in the CR-CRN affected by GoF of KRAS under the action of DBF when also a reaction modeling drug degradation is included in the model. The figure has been designed as panel (A) of Figure 6. A logarithmic scale is used on the time axis; c_D indicated DBF initial concentration.

FIGURE 8

FIGURE 8. Effects of a time-dependent perturbation through the administration of TMT. Similarly to Figure 7, the plot shows the dynamics of the concentration of p-p-ERK in the network affected by a GoF of KRAS where c_D = 50 nM of DBF have been administrated alone (blue line) or together with c_T = 240 nM of TMT given concurrently (green line) or after 12 min, i.e., when the concentration was close to the value at the physiological equilibrium (pink line). The value at the physiological equilibrium is used as reference (black dotted line). Only the degradation of DBF was taken into account.

5 Conclusion and future directions

In this work we have first reviewed the most fundamental aspects of a recently proposed CRN which simulates the behaviour of the signaling network inside a colorectal mutated cell (Sommariva et al., 2021a; Sommariva et al., 2021b). We have described the general context where the CRN is applied, the basic principles underlying the development of a system of ODEs representing the chemical reactions of the network, and we have reviewed the most fundamental properties of the ODEs applied in the simulations.

Next, we have developed new applications concerning comparison between a physiological (healthy) network and a few similar networks resulting from GoF mutations of KRAS, and LoF mutations of PTEN, frequently observed in CRC. Also, drug loaded networks associated with either single or combined targeted drugs have been investigated and compared. The results obtained have been validated using literature data.

The basic novelty of our approach is given by the interaction between global and local aspects in the treatment of mutations and the action of drugs. The most interesting example is concerned with the abnormal value of activated ERK (p-p-ERK) resulting from a GoF of KRAS; unlike other approaches (Pappalardo et al., 2016; Hamis et al., 2021) it is found from the simulations that activated ERK is considerably reduced by administration of DBF, a B-Raf targeting drug. The local time course of the complex p-p-ERK, which is of fundamental interest because of its biologic consequences (Guo et al., 2020; Lavoie et al., 2020; Hamis et al., 2021; Sugiura et al., 2021), is obtained through the use of three different, global, and interconnected networks (physiological, mutated, drug loaded) to assess the general framework which provides both the system of ODEs to be solved, and the required initial conditions.

The methods that have been proposed in this paper may be applied to the prediction of quantitative effects of targeted drugs, and to the optimization of combination therapies for the mutated cell of other cancer types, under rather general conditions.

Our approach focuses on simulations of the signaling network and on the modifications induced by mutations and drugs on the related equilibrium concentrations. For now, we have not considered the cellular response to these changes of the equilibrium state, even though we know that, in general, different equilibrium states correspond to different cell behaviours, as it happens for mutated and physiological equilibrium. In particular, we have not investigated cell behaviours possibly associated with the considered partial GoF of KRAS. Consistently, partial mutations have been defined based on formal assumptions such as a proportional reduction of the values of the rate constants in the physiological network. In our opinion this is only a first, necessary step toward a more comprehensive model of cell behaviour.

In our analysis we have assumed that the parameters of the model are fixed and given, and that a unique equilibrium state exists for every stoichiometric surface. These points need for further investigation. For example, a sensitivity analysis is required to first identify those parameters that are most influent on the equilibrium values, and then to design proper biological experiments to refine their values.

Also, the basic model may require adjustments in order to account for specific effects as the natural process of degradation of drugs and other proteins, or changes of the interactions between proteins, possibly occurring as an answer of the network to modifications induced by drugs. More in general, the model could be extended so as to account for changes induced on the kinetic parameters of the network by the external environment. As an example, the current model mimics the behaviour of a single colorectal cell during the G1-S transition phase. A possible future work may aim at extending this model by considering e.g. selection mechanisms induced on a group of mutated cells by the external environment. This extension may allow our model to explain also more recent experimental evidences by highlighting the importance of the order of mutations on cancer diseases progression (Levine et al., 2019). Furthermore, such an extension may also help to better model the effects of combined therapies. Indeed our simulations seem to suggest that there are no combined concentrations of DBF and TMT that are better than the optimal dose of either alone. This indirectly supports the view that one of the benefits of the combination of two (or more) drugs consists in delaying the onset of mechanisms of drug resistance, that, instead, have been commonly observed in monotherapies (Morkel et al., 2015; Hamis et al., 2021). Future work will be devoted to modeling this type of mechanism.

Finally, we should model glucose metabolism and cell apoptosis. In a sense, our approach to drug action is almost opposite to that of inducing cell apoptosis, in that our main aim has been to restore a (nearly) healthy state, instead of triggering cell death.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: https://github.com/theMIDAgroup/CRC_CRN.

Author contributions

SS, GC, FB, and MP conceived the paper and designed the computational tests. SB and GB realized the simulations. SS, MP, GC, and SB wrote the text of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

SB was granted a scholarship by Roche S.p.A., Italy. SS, MP, and FB have been partially supported by Gruppo Nazionale per il Calcolo Scientifico (GNCS-INdAM). MP acknowledges the support of the “Strengthening a person-centred ecosystem for the co-creation of digital health services for a smart community (DHEAL-COM)” project (PNC E3-2022-23683267).

Acknowledgments

Our colleague Gabriele Zoppoli is kindly acknowledged for discussion on the results of the paper.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsysb.2023.1207898/full#supplementary-material

Footnotes

¹https://www.uniprot.org/uniprotkb/P01116/entry.

²https://www.uniprot.org/uniprotkb/P15056/entry.

References

Anderson, M. W., Moss, J. J., Szalai, R., and Lane, J. D. (2019). Mathematical modeling highlights the complex role of AKT in TRAIL-induced apoptosis of colorectal carcinoma cells. IScience 12, 182–193. doi:10.1016/j.isci.2019.01.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Armaghany, T., Wilson, J. D., Chu, Q., and Mills, G. (2012). Genetic alterations in colorectal cancer. Gastrointest. Cancer Res. 5, 19–27.

PubMed Abstract | Google Scholar

Baharudin, R., Ab Mutalib, N.-S., Othman, S. N., Sagap, I., Rose, I. M., Mohd Mokhtar, N., et al. (2017). Identification of predictive dna methylation biomarkers for chemotherapy response in colorectal cancer. Front. Pharmacol. 8, 47. doi:10.3389/fphar.2017.00047

PubMed Abstract | CrossRef Full Text | Google Scholar

Bedard, P. L., Hyman, D. M., Davids, M. S., and Siu, L. L. (2020). Small molecules, big impact: 20 years of targeted therapy in oncology. Lancet 395, 1078–1088. doi:10.1016/s0140-6736(20)30164-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Berra, S., La Torraca, A., Benvenuto, F., and Sommariva, S. (2022). A fast and convergent combined Newton and gradient descent method for computing steady states of chemical reaction networks. arXiv preprint arXiv:2212.14252

Google Scholar

Biddau, G., Caviglia, G., Piana, M., and Sommariva, S. (2023). Ssi: A statistical sensitivity index for chemical reaction networks in cancer, 2023–2101.bioRxiv. doi:10.1101/2023.01.12.523784

CrossRef Full Text

Biller, L. H., and Schrag, D. (2021). Diagnosis and treatment of metastatic colorectal cancer. Jama 325, 669–685. doi:10.1001/jama.2021.0106

PubMed Abstract | CrossRef Full Text | Google Scholar

Broyde, J., Simpson, D. R., Murray, D., Paull, E. O., Chu, B. W., and Tagore, S. (2021). Oncoprotein-specific molecular interaction maps (sigmaps) for cancer network analyses. Nat. Biotechnol. 39, 215–224. doi:10.1038/s41587-020-0652-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Bruford, E. A., Braschi, B., Denny, P., Jones, T. E., Seal, R. L., and Tweedie, S. (2020). Guidelines for human gene nomenclature. Nat. Genet. 52, 754–758. doi:10.1038/s41588-020-0669-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Conradi, C., and Mincheva, M. (2014). Catalytic constants enable the emergence of bistability in dual phosphorylation. J. R. Soc. Interface. 11, 20140158. doi:10.1098/rsif.2014.0158

PubMed Abstract | CrossRef Full Text | Google Scholar

Castagnino, N., Maffei, M., Tortolina, L., Zoppoli, G., Piras, D., Nencioni, A., et al. (2016). Systems medicine in colorectal cancer: from a mathematical model toward a new type of clinical trial. WIREs Mech. Dis. 8, 314–336. doi:10.1002/wsbm.1342

PubMed Abstract | CrossRef Full Text | Google Scholar

Chellaboina, V., Bhat, S. P., Haddad, W. M., and Bernstein, D. S. (2009). Modeling and analysis of mass-action kinetics. IEEE Control Syst. Mag. 29, 60–78. doi:10.1109/MCS.2009.932926

CrossRef Full Text | Google Scholar

Chen, D., Liu, X., Yang, Y., Yang, H., and Lu, P. (2015). Systematic synergy modeling: understanding drug synergy from a systems biology perspective. BMC Syst. Biol. 9, 56–10. doi:10.1186/s12918-015-0202-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Conradi, C., and Flockerzi, D. (2012). Multistationarity in mass action networks with applications to ERK activation. J. Math. Biol. 65, 107–156. doi:10.1007/s00285-011-0453-1

PubMed Abstract | CrossRef Full Text | Google Scholar

De La Cruz, R., and Kreft, J.-U. (2018). Geometric mean extension for data sets with zeros. arXiv preprint arXiv:1806.06403

Google Scholar

De Martino, A., De Martino, D., Mulet, R., and Pagnani, A. (2014). Identifying all moiety conservation laws in genome-scale metabolic networks. PloS one 9, e100750. doi:10.1371/journal.pone.0100750

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, J. Y., and Richardson, B. C. (2005). The MAPK signalling pathways and colorectal cancer. lancet Oncol. 6, 322–327. doi:10.1016/s1470-2045(05)70168-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Feinberg, M. (1987). Chemical reaction network structure and the stability of complex isothermal reactors-I. The deficiency zero and deficiency one theorems. Chem. Eng. Sci. 42, 2229–2268. doi:10.1016/0009-2509(87)80099-4

CrossRef Full Text | Google Scholar

Guo, Y.-J., Pan, W.-W., Liu, S.-B., Shen, Z.-F., Xu, Y., and Hu, L.-L. (2020). Erk/mapk signalling pathway and tumorigenesis. Exp. Ther. Med. 19, 1997–2007. doi:10.3892/etm.2020.8454

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamis, S. J., Kapelyukh, Y., McLaren, A., Henderson, C. J., Roland Wolf, C., and Chaplain, M. A. (2021). Quantifying erk activity in response to inhibition of the brafv600e-mek-erk cascade using mathematical modelling. Br. J. Cancer 125, 1552–1560. doi:10.1038/s41416-021-01565-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Horst, D., Chen, J., Morikawa, T., Ogino, S., Kirchner, T., and Shivdasani, R. A. (2012). Differential wnt activity in colorectal cancer confers limited tumorigenic potential and is regulated by mapk signaling. Cancer Res. 72, 1547–1556. doi:10.1158/0008-5472.can-11-3222

PubMed Abstract | CrossRef Full Text | Google Scholar

Ingalls, B. P. (2013). Mathematical modeling in systems biology: An introduction. China: MIT press.

Google Scholar

Jin, H., Wang, L., and Bernards, R. (2023). Rational combinations of targeted cancer therapies: background, advances and challenges. Nat. Rev. Drug Discov. 22, 213–234. doi:10.1038/s41573-022-00615-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Jonker, D. J., O'Callaghan, C. J., Karapetis, C. S., Zalcberg, J. R., Tu, D., and Au, H.-J. (2007). Cetuximab for the treatment of colorectal cancer. N. Engl. J. Med. 357, 2040–2048. doi:10.1056/nejmoa071834

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohn, K. W. (1999). Molecular interaction map of the mammalian cell cycle control and dna repair systems. MBoC 10, 2703–2734. doi:10.1091/mbc.10.8.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

Kondratova, M., Sompairac, N., Barillot, E., Zinovyev, A., and Kuperstein, I. (2018). Signalling maps in cancer research: construction and data analysis. Database 2018. doi:10.1093/database/bay036

PubMed Abstract | CrossRef Full Text | Google Scholar

Krogan, N. J., Lippman, S., Agard, D. A., Ashworth, A., and Ideker, T. (2015). The cancer cell map initiative: defining the hallmark networks of cancer. Mol. Cell. 58, 690–698. doi:10.1016/j.molcel.2015.05.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Lavoie, H., Gagnon, J., and Therrien, M. (2020). Erk signalling: a master regulator of cell behaviour, life and fate. Nat. Rev. Mol. Cell. Biol. 21, 607–632. doi:10.1038/s41580-020-0255-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, Y. T., Tan, Y. J., and Oon, C. E. (2018). Molecular targeted therapy: treating cancer with specificity. Eur. J. Pharmacol. 834, 188–196. doi:10.1016/j.ejphar.2018.07.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Levine, A. J., Jenkins, N. A., and Copeland, N. G. (2019). The roles of initiating truncal mutations in human cancers: the order of mutations and tumor cell type matters. Cancer Cell. 35, 10–15. doi:10.1016/j.ccell.2018.11.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Los, M., Roodhart, J. M., and Voest, E. E. (2007). Target practice: lessons from phase iii trials with bevacizumab and vatalanib in the treatment of advanced colorectal cancer. Oncol. 12, 443–450. doi:10.1634/theoncologist.12-4-443

PubMed Abstract | CrossRef Full Text | Google Scholar

Lowenthal, R. M., and Eaton, K. (1996). Toxicity of chemotherapy. Hematology/Oncology Clin. N. Am. 10, 967–990. doi:10.1016/s0889-8588(05)70378-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Morkel, M., Riemer, P., Bläker, H., and Sers, C. (2015). Similar but different: distinct roles for kras and braf oncogenes in colorectal cancer development and therapy resistance. Oncotarget 6–20800. 20785. doi:10.18632/oncotarget.4750

PubMed Abstract | CrossRef Full Text | Google Scholar

Pappalardo, F., Russo, G., Candido, S., Pennisi, M., Cavalieri, S., Motta, S., et al. (2016). Computational modeling of pi3k/akt and mapk signaling pathways in melanoma cancer. PLoS One 11, e0152104. doi:10.1371/journal.pone.0152104

PubMed Abstract | CrossRef Full Text | Google Scholar

Pommier, Y., Sordet, O., Antony, S., Hayward, R. L., and Kohn, K. W. (2004). Apoptosis defects and chemotherapy resistance: molecular interaction maps and networks. Oncogene 23, 2934–2949. doi:10.1038/sj.onc.1207515

PubMed Abstract | CrossRef Full Text | Google Scholar

Porru, M., Pompili, L., Caruso, C., Biroccio, A., and Leonetti, C. (2018). Targeting kras in metastatic colorectal cancer: current strategies and emerging opportunities. J. Exp. Clin. Cancer Res. 37, 57–10. doi:10.1186/s13046-018-0719-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Puszkiel, A., Noé, G., Bellesoeur, A., Kramkimel, N., Paludetto, M.-N., and Thomas-Schoemann, A. (2019). Clinical pharmacokinetics and pharmacodynamics of dabrafenib. Clin. Pharmacokinet. 58, 451–467. doi:10.1007/s40262-018-0703-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocca, A., and Kholodenko, B. N. (2021). Can systems biology advance clinical precision oncology? Cancers 13, 6312. doi:10.3390/cancers13246312

PubMed Abstract | CrossRef Full Text | Google Scholar

Rosen, L. S., Jacobs, I. A., and Burkes, R. L. (2017). Bevacizumab in colorectal cancer: current role in treatment and the potential of biosimilars. Targ. Oncol. 12, 599–610. doi:10.1007/s11523-017-0518-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Saad El Din, K., Loree, J. M., Sayre, E. C., Gill, S., Brown, C. J., and Dau, H. (2020). Trends in the epidemiology of young-onset colorectal cancer: a worldwide systematic review. Bmc Cancer 20, 288–314. doi:10.1186/s12885-020-06766-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Salvatore, L., Calegari, M. A., Loupakis, F., Fassan, M., Di Stefano, B., and Bensi, M. (2019). Pten in colorectal cancer: shedding light on its role as predictor and target. Cancers 11, 1765. doi:10.3390/cancers11111765

PubMed Abstract | CrossRef Full Text | Google Scholar

Santini, C. C., Longden, J., Schoof, E. M., Simpson, C. D., Jeschke, G. R., and Creixell, P. (2019). Global view of the raf-mek-erk module and its immediate downstream effectors. Sci. Rep. 9, 10865. doi:10.1038/s41598-019-47245-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shinar, G., Alon, U., and Feinberg, M. (2009). Sensitivity and robustness in chemical reaction networks. Siam J. Appl. Math. 69, 977–998. doi:10.1137/080719820

CrossRef Full Text | Google Scholar

Sinicrope, F. A. (2022). Increasing incidence of early-onset colorectal cancer. N. Engl. J. Med. 386, 1547–1558. doi:10.1056/nejmra2200869

PubMed Abstract | CrossRef Full Text | Google Scholar

Sommariva, S., Caviglia, G., and Piana, M. (2021a). Gain and loss of function mutations in biological chemical reaction networks: a mathematical model with application to colorectal cancer cells. J. Math. Biol. 82, 55. doi:10.1007/s00285-021-01607-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Sommariva, S., Caviglia, G., Ravera, S., Frassoni, F., Benvenuto, F., and Tortolina, L. (2021b). Computational quantification of global effects induced by mutations and drugs in signaling networks of colorectal cancer cells. Sci. Rep. 11, 19602. doi:10.1038/s41598-021-99073-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Stefani, C., Miricescu, D., Stanescu-Spinu, I.-I., Nica, R. I., Greabu, M., and Totan, A. R. (2021). Growth factors, pi3k/akt/mtor and mapk signaling pathways in colorectal cancer pathogenesis: where are we now? Ijms 22, 10260. doi:10.3390/ijms221910260

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugiura, R., Satoh, R., and Takasaki, T. (2021). Erk: a double-edged sword in cancer. Erk-dependent apoptosis as a potential therapeutic strategy for cancer. Cells 10, 2509. doi:10.3390/cells10102509

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., and Jemal, A. (2021). Global cancer statistics 2020: globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA A Cancer J. Clin. 71, 209–249. doi:10.3322/caac.21660

CrossRef Full Text | Google Scholar

Tariq, K., Tariq, K., Ghias, K., and Ghias, K. (2016). Colorectal cancer carcinogenesis: a review of mechanisms. Cancer Biol. Med. 13, 120–135. doi:10.20892/j.issn.2095-3941.2015.0103

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiwari, A., Saraf, S., Verma, A., Panda, P. K., and Jain, S. K. (2018). Novel targeting approaches and signaling pathways of colorectal cancer: an insight. Wjg 24, 4428–4435. doi:10.3748/wjg.v24.i39.4428

PubMed Abstract | CrossRef Full Text | Google Scholar

Tortolina, L., Duffy, D. J., Maffei, M., Castagnino, N., Carmody, A. M., and Kolch, W. (2015). Advances in dynamic modeling of colorectal cancer signaling-network regions, a path toward targeted therapies. Oncotarget 6, 5041–5058. doi:10.18632/oncotarget.3238

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, C. W.-K., and Lui, R. N. (2022). Early-onset colorectal cancer: current insights and future directions. Wjgo 14, 230–241. doi:10.4251/wjgo.v14.i1.230

PubMed Abstract | CrossRef Full Text | Google Scholar

Xi, Y., and Xu, P. (2021). Global colorectal cancer burden in 2020 and projections to 2040. Transl. Oncol. 14, 101174. doi:10.1016/j.tranon.2021.101174

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y.-H., Chen, Y.-X., and Fang, J.-Y. (2020). Comprehensive review of targeted therapy for colorectal cancer. Sig Transduct. Target Ther. 5, 22. doi:10.1038/s41392-020-0116-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yamamoto, T., Ebisuya, M., Ashida, F., Okamoto, K., Yonehara, S., and Nishida, E. (2006). Continuous ERK activation downregulates antiproliferative genes throughout G1 phase to allow cell-cycle progression. Curr. Biol. 16, 1171–1182. doi:10.1016/j.cub.2006.04.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, P. Y., and Craciun, G. (2018). Mathematical analysis of chemical reaction systems. Isr. J. Chem. 58, 733–741. doi:10.1002/ijch.201800003

CrossRef Full Text | Google Scholar

Zhong, L., Li, Y., Xiong, L., Wang, W., Wu, M., and Yuan, T. (2021). Small molecules in targeted cancer therapy: advances, challenges, and future perspectives. Sig Transduct. Target Ther. 6, 201. doi:10.1038/s41392-021-00572-w

CrossRef Full Text | Google Scholar

Zhu, G., Pei, L., Xia, H., Tang, Q., and Bi, F. (2021). Role of oncogenic kras in the prognosis, diagnosis and treatment of colorectal cancer. Mol. Cancer 20, 143–217. doi:10.1186/s12943-021-01441-4

PubMed Abstract | CrossRef Full Text | Google Scholar