1. The challenges of precision oncology

Over the past decades, two key advances have changed the landscape of cancer research: a predominantly reductionist approach has given way to the widespread use of omics techniques ^[1][2][3], while the therapeutic armamentarium, previously based on chemotherapy and few hormonal drugs, has been enriched with a large variety of targeted drugs (drugs whose molecular targets are precisely defined ^[4]), affecting cancer intracellular signaling or the microenvironment and immune system.

Complexity ^[5] and heterogeneity ^[6] of tumors has prompted the development of precision oncology, aiming to tailor treatments to individual patients ^[7][8] by matching targeted drugs to the molecular alterations that are present in individual tumors. This requires a deep biological characterization of tumors, including omics studies. A fundamental objective is predicting the individual response to treatments, relying on the presence of the target (often a mutant, overexpressed, or fusion protein), but also on predictors formed by panels of biomarkers ^[9]. Another key objective is deciphering the mechanisms of resistance to targeted drugs, such as alterations in the target itself (e.g., secondary mutations), mutations in molecules downstream of the target, and network adaptation mechanisms, including facilitation of resistance by kinase dimerization and feedback mechanisms ^[10][11][12] or activation of parallel pathways ^[13]. A further objective is to identify optimal drug combinations that increase efficacy and overcome resistance ^[14].

The complexity of intracellular signaling and cell-cell interactions makes it particularly challenging to pursue these objectives. Here we review current systems biology methods and provide examples of how they can be valuable tools to tackle this task.

2. Systems Biology

Systems biology studies the collective behavior of different types of molecules involved in a biological process, aiming to reconstruct the system behavior. Systems of many different molecules have behaviors that cannot be simply deduced from the properties of their constitutive elements, and are characterized by “emergent properties” ^{[15][16][17][18]}.

Biological entities are dynamical systems, that evolve in space and time. A snapshot of the system at a certain time point, showing the spatial disposition and the concentrations of molecules, is called a state. The system evolution through different states can be represented graphically as a trajectory in the so-called state space. The aim of systems analysis is to describe the evolution of the system following a perturbation, such as drug treatment.

Dynamical systems evolve towards specific states, called attractors. These may be fixed, steady states, or cyclic trajectories subtending oscillating behaviors, or irregular chaotic trajectories. A system can have more than one attractor, e.g., two different stable steady states (bistability), and evolve toward one or the other depending on the starting state of the system.

Typical emergent properties of dynamical systems are robustness, i.e., the ability to maintain their functions despite perturbations, and adaptation, i.e., the ability to adjust their behavior in response to environmental changes. Robustness and adaptation characterize tumors as complex systems, which almost inevitably adapt to anticancer drugs and develop resistance.

An array of statistical and mathematical modeling techniques can be applied to describe, with different levels of accuracy, dynamical biological processes ^[19].

3. Statistical Methods

Statistical models aim to find associations among variables. Supervised methods deal with different predefined classes of objects (e.g., responders versus non-responders) and try to identify a set of variable values that help us discriminate among classes. Unsupervised methods consider the whole set of data without prior classification and aim to identify intrinsic subsets in the data. Supervised methods include the different types of regression models: linear, logistic, Cox proportional hazards, as well as the so called robust regressions. Unsupervised methods include cluster analysis, aiming to identify relevant subsets in a dataset, and principal component analysis, pointing to reduce the number of variables by combining original variables into a few new, condensed variables. Partial least squares regression makes such variable reduction in a supervised context ^[20].

Machine learning is the part of artificial intelligence (AI) that uses computer algorithms to analyze big datasets to generate predictive models. These algorithms employ statistical tools, both supervised and unsupervised, and are capable to iteratively self-adjust to optimize the performance ^[21].

4. Network Analysis

Another tool to decipher high throughput data is reconstructing networks of interconnected molecules, such as gene and protein networks ^[22] (Figure 1). Network analysis aims first to define the structure or topology of a network, which is linked to its functional properties.

Figure 1. Interaction network of MAPK1. The network is built using the STRING database ^[23]. The top 20 first neighbours and top 20 second neighbours of MAPK1 (aka ERK1) are shown.

Biological networks can be reconstructed following a "bottom-up" (or knowledge-based) approach, whereby the selection of molecules and interactions to be included in the network is based on information extracted from the literature and public databases. Alternatively, a top-down (or data-driven) approach can be followed, reconstructing the network directly from experimental data, for example from omics studies, a process called reverse engineering ^[24]. Different tools are used to reconstruct the interactions among molecules, including statistical tools in Bayesian networks, tools from information theory, or dynamical system-based approaches such as modular response analysis ^[25].

5. Logic Models

Logic models represent the interaction between two molecules as a logical statement ^[26]. In Boolean logic, each molecular species can be in one of two possible states (true or false, on or off, 1 or 0). A Boolean network is formed by a set of Boolean variables, where the connections among variables are defined by Boolean functions. The latter are represented by couples of attributes, specifying the states of the two variables. An activating signal may be represented as 1/1, meaning that when the first node is activated (1, on), also the second becomes activated. An inhibitory signal may be represented as 1/0, meaning that when the first node is activated, the second is inhibited (0, off). Logic modeling does not require knowledge of the detailed mechanistic relationship between nodes, but simply represents the direction and the type (e.g., activating or inhibitory) of relationship. Although logic models mainly yield qualitative results, they can reproduce the evolution of a system. They are suitable to predict the effects of perturbations, such as mutations or exposure to drugs, on a system’s behavior.

6. Mechanistic Models

Mechanistic models that are based on systems of differential equations require a more detailed knowledge of the biological process studied. Building such models requires (i) establishing which molecules and interactions to consider, (ii) choosing mathematical equations to describe each interaction, (iii) finding suitable values for the parameters involved, and then (iv) solving the equations to simulate the behavior of the system, making predictions about system responses to perturbations ^[16][26].

The simplest form of differential equations are ordinary differential equations (ODEs), which use the mean field approximation and are suitable when the number of molecules of each species is large. Modeling biological systems with ODEs involves considering all the processes that modify the level or activity of each relevant molecule. These processes are characterized by rate parameters, such as kinetic constants, whose values must be measured or estimated experimentally. Solutions to ODEs are functions of time, represented as plots on Cartesian axes depicting the evolution of the state of the system over time.

If the assumption of high molecule numbers (more than 1000) is not fulfilled, stochastic differential equations must be used, which consider a probability distribution, instead of a numeric value, for each variable ^[27].

A precise mechanistic modeling requires to represent each different state of any single molecule as an individual variable. Often protein molecules have multiple interacting domains and phosphorylation states, which leads to an exponential increase in the number of variables and ODEs. To overcome this issue, the approach of rule-based modeling has been developed ^[28].

7. Emerging Network Properties Captured by Differential Equation Models

Signal transduction in most cellular networks occurs via cascades of protein phosphorylation-dephosphorylation cycles (Figure 2A), controlled by opposing enzymes, a protein kinase and a phosphatase (Figure 2B,C). These cycles show peculiar behaviors like ultrasensitivity, whereby, when the converting enzymes operate near saturation, the response to an input becomes abrupt ^[29].

Figure 2. A cascade of protein phosphorylation-dephosphorylation cycles.

A cascade of protein phosphorylation-dephosphorylation cycles.

A fundamental feature of networks are feedback loops. Positive feedbacks amplify the signal, whereas negative feedbacks attenuate it. But feedbacks can also favor the occurrence of instabilities, leading to a radical change in the state of the system, known as bifurcations ^[14]. Too strong, long negative feedbacks induce damped or sustained oscillations ^[30] (Figure 3A). Positive feedbacks can cause bistability and hysteresis, where the threshold for jumping from one steady state to the other differs depending on the direction of change of an external signal or parameter (Figure 3B). Hysteresis prevents the easy reversal of a system state, committing to its sustainability. Positive feedbacks occurring in combination with negative feedbacks can give rise to sustained oscillations, called relaxation oscillations, typically observed in cell cycle regulation ^[31].

Figure 3. Oscillations and hysteresis brought about by negative and positive feedback loops.

Under some conditions, hysteresis gives way to an irreversible switch to one of the two steady states ^[32]. Hysteresis and irreversible switches allow the control of multiple irreversible transitions in cellular processes, such as those occurring at cell cycle checkpoints.

8. Modeling Spatiotemporal Network Behavior by Partial Differential Equations

The use of ODE models requires the fulfillment of the assumptions that the molecules are evenly distributed in the modeled compartment, and this condition is not fulfilled for molecules bound to membranes or molecular scaffolds. Often, opposing enzymes, such as a kinase and phosphatase in signal transduction cascades, are localized to different cellular compartments (e.g. a kinase resides at the cell membrane, whereas a phosphatase is homogeneously distributed in the cytosol). This results in the emergence of gradients of phosphorylated and unphosphorylated forms of the protein substrate. Because in this and similar cases, the variables (concentrations) depend not only on time but also on the spatial coordinates, describing the rate of change of these concentrations requires the so called “partial differential equations” (PDEs), which model the spatiotemporal behavior of molecular species. Recently these have been used to model how the biochemical machinery of RhoA GTPase govern cell movements, faithfully reproducing the cellular behavior observed in in vitro experiments ^[33].

9. Mechanistic Models Help Us Understand Resistance to Targeted Therapies

Cellular network adaptations are involved in the development of resistance to targeted therapies. A classic example is resistance to BRAF inhibitors in patients with melanoma. BRAF inhibitors induce conformational changes in BRAF, favoring the formation of RAF homo- and hetero-dimers. If only one RAF protomer is inhibited in a RAF-dimer, the inhibited protomer allosterically activates the other RAF protomer, leading to paradoxical activation of the downstream MAPK pathway. Dynamical models suggest that upon RAF dimerization, for thermodynamic reasons, the affinity of a BRAF inhibitor for one of the RAF protomers increases, but the affinity for the other protomer sharply decreases. This favors the formation of RAF dimers where only one protomer is bound to the inhibitor, with consequent allosteric activation of the inhibitor-free protomer and sustained MAPK signaling activity ^[11][34]. These models predict that a combination of two structurally different inhibitors, binding to different RAF conformations, can overcome resistance, and allow to estimate the levels of synergy or antagonism between RAF inhibitors and MEK inhibitors. These predictions were confirmed with experiments in tumor cell lines ^[12].

Negative feedbacks have been held responsible for the development of acquired resistance to many targeted drugs. Yet, a systematic analysis of network adaptation mechanisms ^[12] shows that, following drug inhibition, negative feedbacks can lead to a transient reactivation or overshoot but cannot fully restore output signaling (Figure 4A,B). Rather, this work highlights three ways in which treatment with a targeted inhibitor can be followed by a complete recovery of pathway signaling: (i) when there are at least two connection routes, activating and inhibitory, from an inhibited upstream drug target to the downstream output, creating a feedforward loop in addition to the other connection route; (ii) when there is a crosstalk among two pathways, whereby inhibition of one of them can favor the activation of the other; (iii) when a kinase inhibitor induces kinase dimerization, as described above for RAF inhibitors ^[12]. Combined with kinase dimerization, negative feedbacks facilitate drug resistances.

10. Signaling Network Models Can Predict Drug Sensitivity

Modeling the dynamics of signaling pathways to drug perturbations can predict cell type–specific responses to small-molecule therapeutics and prioritize primary drug targets and their combinations. Dynamic logic models of signaling networks were built for colorectal cancer cell lines based on a large-scale signaling perturbation screening ^[35]. Simulated signaling dynamics correlated with some drug sensitivities, and a drug combination predicted to overcome resistance to MEK inhibitors by co-blockade of GSK3 was validated experimentally. A dynamic model of the estrogen receptor-induced proliferation of MCF-7 breast cancer cells was built to predict responses to endocrine therapy ^[36]. This ODE model was able to predict the responses to a combination treatment of estrogen deprivation and the drug palbociclib inhibiting Cdk4/6. Immune checkpoint inhibitors (ICIs) greatly enhanced cancer treatment, yet many patients are intrinsically resistant to anti-programmed cell death protein 1 (anti-PD1) and anti-cytotoxic T-lymphocyte-associated protein 4 (anti-CTLA4) therapies or they become resistant after initial response. A logic network model encompassing PD1, CTLA4 and other inhibiting and activating checkpoints ^[37] was able to recapitulate results of existing experimental studies of anti-PD1 and anti-CTLA4 therapies, and suggested the best combinations of ICIs with drugs targeting other checkpoints.

11. Patient-Specific Network Modeling

Models can be used to carry out patient-specific network simulations and to construct patient-specific dynamic biomarkers. For this purpose, variables such as the amounts, mutational status and copy number variation of relevant molecular species are measured directly on a patient tumor sample and the values are entered into the model, to perform patient-specific simulations of the system behavior.

Several studies have exploited systems-level models of apoptosis to predict the patient responses to chemotherapy. Patient-specific ODEs-based models of the intrinsic apoptotic pathway, based on the concentrations of five key molecules determined in samples of colorectal tumors, were the only significant predictor of patient outcome at multivariate analysis and outperformed predictors based on statistical analyses of apoptotic molecules ^[38]. A logic model of the intrinsic and extrinsic apoptosis pathways was used to simulate the apoptotic response to different drugs or drug combinations in pancreatic tumor samples and cell lines, allowing to predict effective drug combinations that were confirmed experimentally ^[39].

Other studies have been performed to define the best combination therapy according to the set of driver genes of an individual tumor ^[40], to reconstruct metabolic networks ^[41] and to investigate the metabolic differences subtending different tumor phenotypes, such as resistance or sensitivity to radiation therapy ^[42].

The c-Jun N-terminal kinase (JNK) pathway is a MAPK cascade mediating apoptosis in response to different types of stress, including chemotherapeutic agents. JNK may undergo either a gradual activation in response to growth factors, promoting cell survival and proliferation, or an ultrasensitive, switch-like activation in response to stress, leading to apoptosis. Reconstruction of the JNK pathway in neuroblastoma cell lines identified a positive feedback from JNK to MKK7 as responsible for the ultrasensitive switch-like apoptotic response ^[43]. Patient-specific simulations obtained by an ODEs-based model of JNK pathway showed that JNK response is increasingly impaired with increasing stage of the disease and has independent prognostic value in multivariate analysis ^[43].

12. Conclusions

Despite major advances, cancer treatment remains an open challenge in many respects. Identifying the driver molecular alterations in a tumor is only part of the solution ^[44], and to obtain long-term control, or even cure, of complex tumors will likely require that (i) their key molecular alterations are correctly targeted, and (ii) network adaptations causing drug resistance are prevented by properly designated drug combinations. Describing and analyzing the intrinsic behavior of the biological processes that underlie tumor pathology, dynamical models will potentially allow to identify the key points for effective interventions with target drugs, substantially delaying or preventing resistance.