Optimizing Dosage-Specific Treatments in a Multi-Scale Model of a Tumor Growth

The emergence of cell resistance in cancer treatment is a complex phenomenon that emerges from the interplay of processes that occur at different scales. For instance, molecular mechanisms and population-level dynamics such as competition and cell–cell variability have been described as playing a key role in the emergence and evolution of cell resistances. Multi-scale models are a useful tool for studying biology at very different times and spatial scales, as they can integrate different processes occurring at the molecular, cellular, and intercellular levels. In the present work, we use an extended hybrid multi-scale model of 3T3 fibroblast spheroid to perform a deep exploration of the parameter space of effective treatment strategies based on TNF pulses. To explore the parameter space of effective treatments in different scenarios and conditions, we have developed an HPC-optimized model exploration workflow based on EMEWS. We first studied the effect of the cells’ spatial distribution in the values of the treatment parameters by optimizing the supply strategies in 2D monolayers and 3D spheroids of different sizes. We later study the robustness of the effective treatments when heterogeneous populations of cells are considered. We found that our model exploration workflow can find effective treatments in all the studied conditions. Our results show that cells’ spatial geometry and population variability should be considered when optimizing treatment strategies in order to find robust parameter sets.

1 Introduction

Optimizing drug treatment and efficiently screening the effect of drugs is key to improving clinical treatments and ultimately extending patients’ life expectancy (Kessler et al., 2014). The emergence of resistant cancer cells is a complex phenomenon due to the inherent complexity of biological (Shaffer et al., 2017), the interplay of processes that occur at different scales, and an environment with an active role in this resistance (Lee et al., 2012; Goldman et al., 2015). Molecular mechanisms and population-level dynamics such as competition and cell–cell variability have been described as playing a key role in the emergence and evolution of cell resistances (Kim et al., 2018). For instance, high gene expression variability has been linked to aggressiveness in chronic lymphocytic leukemia (Ecker et al., 2015). Genetic heterogeneity and phenotype variability have also been related to the emergence of cell resistance (McGranahan and Swanton, 2015; Brady et al., 2017; Shaffer et al., 2017). Furthermore, the environment has been described to have an effect on the cells’ response to drugs: 2D-cultured cell line screens failed in clinical studies (Horvath et al., 2016) as cell cultures do rarely recapitulate the heterogeneity and drug sensitivity of the original tumor (Jabs et al., 2017).

Multi-scale models (MSM) are a useful tool for studying biology at very different time (no s) and spatial scales, as they can integrate different processes occurring at the molecular, cellular, and intercellular levels (Metzcar et al., 2019; Montagud et al., 2021). In the domain of cancer biology, MSMs have been used to connect cellular mechanisms underlying cancer drug resistance to population-level patient survival (Sun et al., 2016), study the role of physiologic resistance due to diffusion gradients of different nutrients and drugs (Frieboes et al., 2009), and quantitatively characterize pressure for invasion (Anderson et al., 2006), among many other applications (Metzcar et al., 2019). In general, multi-scale models provide a genotype-to-phenotype simulation framework, which is ideal for the study of in silico drug screenings (Flobak et al., 2015), the optimization of treatment regimens (Akasiadis et al., 2021), and the exploration of genetic or environmental perturbations (Letort et al., 2018).

Multi-scale simulation can be used to conduct in silico experiments and generate new experimentally testable hypotheses, accelerating the discovery of new potential treatment strategies (An, 2010). Nevertheless, due to the hybrid approaches used to describe multi-scale models (e.g., discrete, continuous, and stochastic), these models cannot be studied using formal analytical tools, and thus the analysis and exploration of simulated trajectories require complex workflows to guide the exploration of the parameter spaces associate with these models (Ozik et al., 2016). For this reason, distributed workflows to perform parallel optimization via simulation and model exploration are critical tools for exploiting the full potential of simulations (Ozik et al., 2018b; Reuillon et al., 2013). Model exploration workflows are required to efficiently fit parameters for which there are no available experimental measurements (Akasiadis et al., 2021; Ozik et al., 2019), explore complex and vast parameter spaces, and optimize user desirable goals, such as the space of optimal treatment strategies for a given cancer model. Optimization methods such as evolutionary algorithms have proven their usefulness in such studies for fitting unknown parameters (Akasiadis et al., 2021), as well as high-throughput hypotheses testing in cancer research (Ozik et al., 2018a).

In previous work, Letort et al. (2018) developed the multi-scale model of 3T3 fibroblast spheroids that integrates the Cell Fate Boolean network Calzone et al. (2010) inside individual cell agents. The Boolean network rules the phenotype of the cells (e.g., proliferation and apoptosis) based on the environmental conditions (e.g., drugs presence and oxygen concentration). The authors used the model to investigate the tumor response to different regimes with tumor necrosis factor supplies (TNF) and reported complex behaviors in the simulated conditions. While a set of values of pulse period, pulse duration, and TNF concentration was optimal to reduce the number of alive tumor cells, different sets of values turned the cells resistant to TNF (Letort et al., 2018). The effects of TNF in the Boolean model reported by Calzone et al. (2010) are multifaceted: TNF triggers cells to go from a naive to a proliferative state and commits cells to necrosis and apoptosis. Once the cells are committed to either survival, necrosis, or apoptosis, they cannot go back, causing resistance due to phenotypic variability. Interestingly, it has been described that prolonged TNF exposure causes the cells to be resistant to the effect of the cytokine (Lee et al., 2016).

In the present work, we use an extended hybrid multi-scale model to perform a deep exploration of the parameter space of effective treatment strategies based on TNF pulses to unravel the mechanistic details behind the complex emergent dynamics of the TNF pulses in in silico experiments and guide the optimization of effective treatments. We extended the multi-scale model of 3T3 fibroblast spheroid by integrating an explicit kinetic description of the TNF-receptor dynamics based on the molecular biology of the TNF receptor (Fischer et al., 2011; Li et al., 2013; Sedger and McDermott, 2014). Furthermore, we couple the TNF-receptor kinetic model with the cancer cell Boolean model from Calzone et al. (2010) to simulate the downstream propagation of the signal that induced the binding of the TNF. To explore the parameter space of effective treatments in different scenarios and conditions, we have developed an HPC-optimized model exploration workflow based on EMEWS (Ozik et al., 2018b). Our workflow includes two previously used model exploration strategies, sweep search and genetic algorithm (Akasiadis et al., 2021; Ozik et al., 2019), together with a new approach named the Covariance Matrix Adaptation Evolutionary Strategy, which exhibited good convergence in global optimization problems with continuous variables (Hansen and Ostermeier, 2001).

We applied our framework to characterize the space of effective treatments in different experimental scenarios by simulating the treatment outcome with our multi-scale model of tumor growth. We first studied the effect of the cells’ spatial distribution in the values of the treatment parameters by optimizing the supply strategies in 2D monolayers and 3D spheroids of different sizes. We found that our model exploration workflow can find non-trivial in silico drug scheduling strategies that minimize the tumor below 1% of its initial size while avoiding the emergence of resistant cells. Our results also show that effective treatment strategies can be found in the two different cell geometries studied. We also found that the parameter spaces of effective treatments for the 2D monolayer and 3D spheroid exhibit different distributions for the parameters. We later study the robustness of the effective treatments when heterogeneous populations of cells are considered. Specifically, we model population heterogeneity by introducing different levels of cell-based variability into the kinetic parameters of the TNF-receptor models. The parameters’ variability aims to mimic population-level variability in the kinetic parameters of the receptor, as well as different levels of expression in the receptor, among different cells. We found that effective treatment strategies are robust to a low level of variability, whereas, with a high level of variability, those treatment strategies optimized for populations with no variability cannot reduce tumor growth. However, when the treatments are optimized directly on a heterogeneous population, we observe that the optimization algorithms can retrieve effective treatment.

Altogether, we found that our model exploration workflow can find effective treatments in all the studied conditions, showing that multi-scale simulations and model exploration are promising tools for in silico exploring treatment strategies. Finally, our results also show that cells’ spatial geometry and population variability should be considered when optimizing treatment strategies to find robust parameter sets. In future work, we plan to extend these results by studying other experimental setups and different cancer models.

2 Materials and Methods

2.1 Hybrid Multi-Scale Model of Cancer Cells With Signaling

Herein, we present a multi-scale model of tumor growth that considers, at the individual cell level, the dynamics of the tumor necrosis factor (TNF) receptor and its downstream effect using a hybrid approach (Figure 1). Our model was implemented using PhysiCell (Ghaffarizadeh et al., 2018) together with the PhysiBoSS add-on (Ponce-de-Leon et al., in preparation). The microenvironment is simulated in both the 2D and 3D domains, and it accounts for the presence of oxygen and the cytokine tumor necrosis factor (TNF). On the contrary, cells are simulated as individual agents, including intracellular submodels that account for the cell cycle, the different death models (i.e., necrosis and apoptosis), a model for TNF receptor dynamics, and a gene regulatory network.

FIGURE 1

FIGURE 1. Diagram representing the intracellular submodels of the multi-scale model of tumor growth. Each individual cell agent has a kinetic model of the TNF receptor dynamics connected to the microenvironment through the presence of surrounding TNF and coupled to the Boolean network through a transfer function. The Boolean network has three readout nodes (proliferation, NonACD, and apoptosis), which rule the fate of the cell agent.

For the cell cycle, we use PhysiCell live cell cycle with a doubling time of 22h, and for the death models, we used PhysiCell standard ones with default parameters. The binding of the TNF to its receptor is modeled using mass-action kinetics in which TNF binds to a cell receptor TNFR at a given rate k_bind; the complex TNF-TNFR is internalized at a rate of k_endo where the TNF is degraded and the receptor recycled at a rate of k_recycle (see Supplementary Figure S1). The TNF-receptor submodel was developed based on the known molecular biology of the molecular system (Fischer et al., 2011; Li et al., 2013; Sedger and McDermott, 2014). The equation below describes the submodel for the TNF-receptor dynamics:

$\begin{array}{cc}\hfill & \frac{\left[R_e\right]}{dt}=-k_{bind}\left[R_e\right]\left[TNF\right]+k_{recycle}\left[R_i^{\ast }\right]\hfill \\ \hfill & \frac{\left[R_e^{\ast }\right]}{dt}=k_{bind}\left[R_i\right]\left[TNF\right]-k_{endo}\left[R_e^{\ast }\right]\hfill \\ \hfill & \frac{\left[R_i^{\ast }\right]}{dt}=k_{endo}\left[R_e^{\ast }\right]-k_{recycle}\left[R_i^{\ast }\right]\hfill \end{array}(1)$

where [R], [TNF], and [R*] are the concentrations of the receptor, TNF, and TNF-TNFR complex, respectively. Furthermore, the TNF-TNFR complex [R*] can be found in two states, in the cell membrane $(R_e^{*})$ or internalized $(R_i^{*})$.

The gene regulatory model used is an extended version of the Boolean network (BN) reported in Calzone et al. (2010) and is simulated using the MaBoSS algorithm. The BN is coupled to the agent in two different ways (see Supplementary Figure S2). The BN has an input node that represents the presence of TNF and is coupled to the amount of active TNF-TNFR complex $(R_e^{*})$ through a transfer function that converts the continuous value of $R_e^{*}$ into a Boolean one. Additionally, the BN has three mutually exclusive output nodes representing three alternative cell fates: proliferation, apoptosis, and NonAD (non-apoptotic death or necrosis). The fate or phenotype of each cell agent is ruled by the current state of the fate nodes of its internal regulatory network. For instance, if the proliferation node is active, the cell will grow and divide, whereas if the apoptosis becomes active, the cell agent will commit to apoptosis (see Supplementary Figure S2).

2.2 Model Exploration Framework

2.2.1 Workflow Overview

The parallel simulation framework used in our evaluation is a workflow that follows the Extreme-scale Model Exploration with Swift (EMEWS) paradigm. It uses the spheroid_TNF_v2 as an example model and is publicly available in our online repository¹. An overview of our model exploration workflow is shown in Figure 2. We have integrated three different search strategies: 1) a sweep search approach that evaluates a predetermined set of candidate parameters (generated from uniform sampling, or a regular grid), 2) a Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES), and 3) a genetic algorithm (GA). For the cases of CMA-ES and GA, we use the available implementations provided by the DEAP package (version 1.3.1), a mature and widely used package for evolutionary optimization (Fortin et al., 2012). Using EMEWS queues, multi-scale simulation instances are configured with the specific parameter values corresponding to the points that each exploration procedure targets and are then submitted for parallel execution in an HPC environment.

FIGURE 2

FIGURE 2. Workflow overview. The diagram depicts the structure of the model exploration workflow. EMEWS communicates to the different search strategies using a queue system. The search strategy generates candidate parameters, and the treatments to be evaluated via PhysiBoSS simulations are distributed as parallel jobs to the HPC infrastructure. Upon completion, the simulation outcomes are returned to EMEWS, which, in the case of the GA and CMA-ES, sends the fitness of the evaluated parameters, so the algorithm can update its internal state and generate new candidate parameters.

The number of each “batch” of points is relative to the number of computational nodes available in the HPC. The multi-scale simulator incorporates PhysiCell (v1.7) (Ghaffarizadeh et al., 2018) together with the PhysiBoSS 2.0 extension (Ponce-de-Leon et al., in preparation), which is an add-on version of PhysiBoSS. We merged our mass-action kinetics model, as explained in Section 2.1, with the multi-scale model proposed by Letort et al. (2018) that used the Boolean model from Calzone et al. (2010). Finally, the simulation results are returned, the points are evaluated according to the performance of the particular drug treatment, and the workflow iterates over the next “batch” of points. Each point is a three-dimensional vector that configures the following simulation parameters: 1) the duration of the TNF pulse; 2) the TNF pulse period; and 3) the concentration of TNF. The exploration space ranges from 5 to 800 min for the pulse period, from 5 to 200 min for the pulse duration, and from 0.001 to 1 ng/L for the TNF concentration.

We checked the number of alive tumor cells at the last time point of each simulation to evaluate the results of each particular treatment and used these values as the fitness or objective of the optimization algorithms. Note that, to ensure that the characterization is robust and not a subject of extreme randomness accruing from inherent PhysiBoSS stochasticity, we perform three replicate simulations with the same configuration parameters, using a different seed to initialize the random number generator, and calculate the average value of the final alive tumor cells count over the replicates as the final score. We now proceed to describe the different search methods we use in detail.

2.2.2 Sweep Search

The sweep search comprises a simple exhaustive approach that requires the user to specify a predetermined number of points to be evaluated. Our code offers a points generating script, which can be configured to choose among different distributions. In other words, points can either be selected to belong on a grid, with equal distances between each point along the dimensions, or a second option is to select random points by sampling particular probability distributions. For the purposes of the experimentation presented in this study, we have implemented the uniform distribution point selection, though this can be easily configured to use other types, such as Gaussian and Beta.

2.2.3 Covariance Matrix Adaptation Evolutionary Strategy

Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES) is a stochastic, derivative-free method for numerical optimization for black-box optimization functions (Hansen and Ostermeier, 2001). This method requires, as input, a set of points, σ value that controls the range of exploration, covariance matrix C used to guide the search, the number of points population to execute the algorithm upon, and a total number of iteration or stop criteria. In a nutshell, CMA-ES generates an initial population sampling from a multivariate normal distribution, evaluates each generated point, and then calculates mutation steps of the best points to form the mutation distribution. By this, a new population of points is generated and evaluated, and the iterative process continues up to a user-defined number of times. Note that, for every update of the mutation distribution in each algorithm iteration, all past paths from previous iterations are also considered, and the most favorable points are granted a larger probability for being selected by the evolutionary strategy. This way, the length of each mutation step can be adapted to be longer in cases of greater fitness score improvement or shorter for the opposite case.

2.2.4 Genetic Algorithm

Genetic algorithm (GA) is a widely known and tested metaheuristic approach, also belonging to the family of evolutionary strategies algorithms, which mimics the evolution principles of biological organisms and operates directly on the values of points (Holland, 1975; Whitley, 1994). Similar to CMA-ES, an initial population of points is generated and evaluated, and then, following an iterative approach, a series of genetic operators are applied to each of them in order to produce the next evolved population of points. More specifically, the GA applies the selection, crossover, and mutation operators. Typically, the first operator selects the evaluated individuals in a weighted manner so that the ones with better fitness scores have an increased probability of being selected to proceed to the next generation, compared to fewer fit points. Then, the crossover mixes the point values in a principle similar to that of the gene propagation from parents to offspring as it happens in organisms. Finally, the mutation operator changes a point value (e.g., one of its dimensions) with a small probability, similarly to the process that has been observed in DNA sequences. The main idea is that combinations of points with good fitness scores would lead to even better ones, especially if the search domain is smooth. However, because the algorithm considers only the previously observed fitness scores, without having any other domain-specific knowledge, as a consequence, its search may be constrained around locally optimal points, never managing to reach the global optimal ones. Despite these shortcomings, GAs have been shown to work very well for non-smooth search spaces (Fitzpatrick and Grefenstette, 1988; Tang et al., 1996).

3 Results

3.1 Multi-Scale Simulations and Model Exploration Setup

Herein, we use a multi-scale model of tumor growth to investigate different treatment strategies. The model, which is also used in Akasiadis et al. (2021), simulates the dynamics of a population of cancer cells growing under different drug treatment conditions. A treatment strategy consists of the supply of periodic pulses of the cytokine tumor necrosis factor (TNF) with fixed duration and concentration (see Section 2). At the molecular level, when the TNF binds to the cell’s receptor TNFR forming a complex, and the TNF-TNFR complex concentration reaches a given threshold, the signal is propagated through the Boolean regulatory network, inducing cell death. However, if the stimulus is sustained for a longer period of time, cells activate the NFkB node and the survival node, becoming resistant to the death induced by the TNF. For this reason, optimal treatments should expose the cell for a sufficient time to induce death but not too much as to become resistant to it (Letort et al., 2018).

To explore the parameter space associated with the treatment, we have extended our model exploration workflow based on EMEWS (Akasiadis et al., 2021). In each in silico experiment, we simulate the growth of a population of cancer cells for 4,640 min (i.e., three days) subject to a given treatment strategy. In order to account for the inherent stochasticity of the model, each simulation is always run in three replicates and the average behavior is considered (see Section 2). We evaluate the effect of the treatment strategies by analyzing the total number of alive cells at the end of the simulations relative to the initial population size and use these values as the score or cost function associated with a treatment strategy. We define, as effective treatments, those strategies that reduce the number of the alive cell below 1% of its initial numbers in the three replicates. Based on this definition, we investigate the parameter space of the effective treatments in two different spatial arrangements of cells: a monolayer disc of radius 100 μm (151 cells) and a 3D spheroid of radius 100 μm (1,173 cells).

3.2 Effective Treatment Parameters Differs for 2D and 3D Cell Arrangements

To investigate the structure of the parameter space of the effective treatments, we perform a uniform sampling of 10,000 candidate sets of parameters corresponding to different treatment strategies. We use these sets of parameters (sweep search) as inputs for multi-scale simulations and evaluate each treatment effect on the growth of the cancer cells in the 2D and 3D arrangements. From the 10,000 evaluated parameter sets in each condition, the results show that 113 strategies are effective treatments for the 2D setup, whereas, in the 3D case, only 11 strategies are effective treatments (see Supplementary Table S1). This indicates that the region containing the effective treatments for the 3D spheroid is more constrained than in the 2D disc arrangement.

Interestingly, if we restrict the definition of effective treatments to kill all cancer cells at the end of the simulation and in the three replicates, only eight sets of parameters can reach the goal for the 2D cases, whereas no effective parameter sets are found for the 3D setup (see Supplementary Table S1 and Supplementary Figure S3).

We compared the distributions and summary statistics of the parameters of the effective treatments in the two arrangements (Figure 3). The comparison shows that the distribution values for the evaluated parameters indicate that the effective treatments of the 3D arrangement are notably more constrained than those that work in the 2D arrangement, in particular regarding the concentration of TNF and the pulse duration. In general, the effective treatments in 2D arrangements exhibit bigger values and larger ranges in the three parameters (Table 1).

FIGURE 3

FIGURE 3. Effective TNF treatment parameters distribution from uniform random sampling. The distribution of the three parameters of the pulse treatment, sampled from a uniform distribution and filtered to belong to the feasible region that can reduce the tumor size below 1% of its initial size.

TABLE 1

TABLE 1. Summary statistics for the parameters from the effective treatments found by sampling 10,000 random candidates.

We also found that some parameters’ combinations exhibit correlations (see Supplementary Figure S4). On the one hand, we observe that the Pulse period positively correlates with TNF in the 2D and 3D but only correlates with pulse duration in the 2D case. On the other hand, TNF shows a negative correlation with pulse duration only for the 2D case. Altogether, these correlations indicate that the treatment parameters can compensate each other: a shorter pulse might be as effective as a longer one if it carries more TNF.

Although herein we are considering a spheroid composed only of tumor cells, more complex scenarios will also contain healthy cells. In such cases, effective treatments will also need to consider the cytotoxic effect of the drug on the healthy cells. To address this issue, we have also calculated the total concentration of TNF supplied during each treatment based on the pulse parameters as follows:

TNF_total = Pulse_duration*Pulse_{concentration}*n,

where n is the total number of pulses supplied calculated by dividing the total treatment (simulation) duration (min) over the pulse period (min). Using the calculated values, we ranked the feasible solutions to find the ones that minimize the total concentration of TNF used during the whole treatment (see Supplementary Table S1). The results show that the distribution of total TNF is biased to lower values of the total TNF supplied in both the 2D and 3D (see Supplementary Figure S5). Strikingly, when we analyze the effective treatments that minimize the total amount of TNF supplied in the 2D and 3D, we found a very similar value of 5 ng/L.

3.3 Optimal Treatment Parameters Differs for 2D and 3D Cell Arrangements

To further investigate the structure of the parameter space of the effective treatments, we conducted an optimization via simulation to find the set of treatment parameters that minimizes the number of alive cells at the end of the application of the treatment. We performed the parameter optimization in both cell arrangements (i.e., 2D disc and 3D spheroid of radius 100 μm), focusing on the same parameters as in Section 3.2: pulse duration, pulse period, and TNF. These optimizations were run using two evolutionary algorithms: Covariance Matrix Adaptation Evolutionary Strategy (CMA-ES) and genetic algorithm (GA) (see Section 2 for details).

The results showed that both algorithms converge to optimal (leaving no remaining cells) or near-optimal sets of parameters for the 2D and 3D case, respectively (Figure 4). Interestingly, in both cases, the GA found effective treatments after two iterations, whereas the CMA-ES only does after 15 iterations. Nonetheless, both algorithms are capable of finding effective treatments. In addition, the CMA-ES algorithm also estimated a multivariate normal distribution for the region of effective treatments parameters that is updated at each iteration. In the last iteration, the population sampled by the CMA-ES showed a very low variance (Figure 4) in the 2D and the 3D arrangements, indicating that the estimated distribution captures, at least, part of the structure of the parameters associated with the effective treatments.

FIGURE 4

FIGURE 4. Algorithmic convergence for the optimization of treatment parameters. The average % of alive cells at the end of each iteration step of the CMA-ES and GA optimization algorithms. Panels (A) and (B) show the convergence of the algorithms for simulations considering a population of cancer cells arranged in a 2D disc and a 3D spheroid, respectively.

We compared the parameter sets for the effective treatments predicted by both algorithms and uncovered that each one converges to different regions of the parameter space (see Supplementary Figure S6). The CMA-ES found effective treatments parameter distributions different for the 2D and 3D. In both cases, the parameter ranges were narrower than those found in uniform sampling. Moreover, the CMA-ES converged to distributions for the pulse period different for the 2D and the 3D arrangement. In addition, the parameters corresponding to effective treatments found by the GA were more scattered showing wider ranges of values. In both cases, the pulse duration seemed to be the less critical or constrained parameter.

We analyzed the time course of the total number of alive, apoptotic, and necrotic cells for one of the optimal treatments found. Furthermore, for this simulation, we also checked the internal state of the TNF-receptor model (i.e., the values for R_e, $R_e^{*}$, $R_i^{*}$ variables) at each time step of the simulation and averaged these values over all the cells to investigate the coarse grain dynamics of the submodel. Figure 5A shows how the number of alive cells periodically drops until it reaches zero as a consequence of the TNF pulses. When analyzing the average dynamics of the TNFR receptor model in effective treatments, we found that the TNF pulse needs to trigger the activation of the receptor of 50% of the population for short periods of time to be effective (5b). If the average number of cells that got activated is lower than this threshold, then the rate of cells entering necrosis will be lower than the population growth rate, and therefore the number of alive cancer cells will steadily grow. On the contrary, if the average number of activated cells is above this threshold for extended periods of time, many cells will become resistant to the treatment. We found similar results for the case of the 3D spheroid (data not shown).

FIGURE 5

FIGURE 5. Time course for effective treatment in 2D cell arrangement. Panel (A) shows the time course for the number of alive, apoptotic, and necrotic cells. The light grey curve shows the TNF pulses. Panel (B) shows the average state of the TNFR receptor model across all the alive cells. The horizontal line indicates the threshold at which the signal induced by the binding of the TNF propagates downstream to the Boolean network.

3.4 Robustness Analysis of the Effective Treatments in Heterogeneous Populations

In the previous section, we showed that effective treatments could be found for 2D and 3D cell arrangements, using either the CMA-ES or the GA algorithms. Nonetheless, those treatments were optimized on monoclonal or homogeneous tumors, that is, the population of cells with identical parameters. In this section, we study the robustness of effective treatments by studying heterogeneous populations of cells. We forced the population heterogeneity by introducing variability into the three kinetic parameters of the TNF receptor, that is, the TNFR binding rate, the TNFR endocytosis rate, and the TNFR recycling rate. The variability is applied by considering a normal distribution centered in each parameter’s default value and with a standard deviation, the control parameter that varies from zero (homogeneous population) to one (almost uniformly distributed random parameters). Then, when the population is initialized, the kinetic parameters of each cell are sampled from the corresponding distribution.

To evaluate the robustness of the effective treatments in heterogeneous populations, we considered the top 30 effective treatments parameter sets that had no final tumor cells in any of their replicates for the 2D and for the 3D cell arrangements. Then, for each set of parameters, we run the simulations with different levels of variability from 0 to 1. As expected, we observed that, for low values of the variability control parameter $(<0.2)$, most of the evaluated effective treatments still can reduce the initial tumor size to the 1% of the initial size. However, for higher variability values $(>0.2)$, most of the evaluated treatments could not reduce the tumor size below its initial size (Figure 6). This indicates that, as the level of variability on the TNFR receptor kinetic parameters increases, the effectiveness of the treatments dramatically decreases (Figure 6).

FIGURE 6

FIGURE 6. Evaluation of the top 30 best effective treatments in heterogeneous populations of cancer cells with different degrees of variability in the kinetic parameters of the TNFR. For a given value kinetic parameters’ variability (x-axis), each pair of boxes depicts the distribution % of alive cells at the end of the simulations obtained after evaluating the top 30 best effective treatments founded when zero variability was considered.

Interestingly, the critical value at which most effective treatments were no longer effective is different for the 2D and the 3D cell arrangements. In the case of the 2D disc, the critical value is close to 0.25, whereas, in the 3D spheroid cases, this value is around 0.15 (Figure 6). When the variability value is above this threshold, some of the sampled kinetic parameters make the cell insensitive to the treatment, and thus it can grow even in the presence of TNF. The differences in the critical threshold were possibly due to the different number of initial cells considered in the 2D and the 3D cell arrangements. To assess this, we tested the effect of the radius size of the 2D simulations in this robustness analysis.

We evaluated radius sizes of 50, 275, and 500 μm corresponding to initial tumor sizes of 37, 1,069, and 3,559 cells, respectively. We tested three levels of variability with each radius size in three replicates and averaged their results. As already discussed, we can observe that the more variability, the worse the outcome (see Supplementary Figure S7). Interestingly, we did not see a clear correlation between the radius length and the decrease of the effectiveness of the treatments, with the 50 μm being the one with worse outcomes with the higher variability level.

3.5 Optimization via Simulation Can Find Effective Treatments in Heterogeneous Populations of Cancer Cells

To evaluate the performance of model exploration workflow in more complex scenarios, we investigated the optimization of treatment strategies in tumors with different levels of heterogeneity. This use case was considered as a way to evaluate what would happen when using this methodology in a less ideal situation as it can be the drug screens in cell lines or tumors with heterogeneous non-clonal cells. For this purpose, we conducted the optimization via simulation to find effective treatments in two conditions with different degrees of variability in the kinetic parameters.

At first, we set the variability value on the kinetic parameters of the TNFR receptor model to 0.25 and then run the GA and CMA-ES to find treatments that minimize the total number of alive cells. We found that, with this level of variability, neither of the algorithms could converge (see Supplementary Figure S8). Nonetheless, this did not prevent the algorithms from finding effective treatments for both arrangements: the 2D disc and the 3D spheroid. While, for the 2D disc, several candidate sets of parameters were found, only a few candidate effective treatments were found by the CMA-ES for the 3D spheroid.

Interestingly, for the 2D cell arrangement, the CMA-ES algorithm was able to find two optimal treatment strategies, that is, a set of parameters for the TNF pulse that kills every tumor cell in all three replicates. These two sets of treatment parameters are similar to those effective treatments that worked when variability was not introduced (see the previous section). Figure 7A shows the time course for the effective treatments optimized in the 2D cell arrangement. The plot shows how the number of alive cells steadily decreases until zero (Figure 7B). It also shows the average internal dynamic of the receptor model exhibiting some noisy behavior due to the heterogeneity present in the population.

FIGURE 7

FIGURE 7. Time course for effective treatment in the 2D cell arrangement with variability of 0.25. Panel (A) shows the time course for the number of alive, apoptotic, and necrotic cells. The light grey curve shows the TNF pulses. Panel (B) shows the average state of the TNFR receptor model across all the alive cells. The horizontal line indicates the threshold at which the signal induced by the binding of the TNF propagates downstream to the Boolean network.

We also performed a similar experiment with a higher variability value (0.50) on the kinetic parameters of the TNFR receptor model. We ran the treatment optimization using the GA and CMA-ES and found that, as expected, with this level of variability, the convergence of both algorithms was even worse than for the case of 0.25 (see Supplementary Figure S9). Furthermore, the algorithms could only find effective treatments for the 2D disc arrangement. Nevertheless, the number of different parameter sets found were fewer than in the previous scenario, with a variability value of 0.25, as expected for a more complex landscape. For the 3D spheroid case, the best treatment reduced the initial tumor size to 1.05% of its initial size when averaged over the three replicates. If we relax the definition of effective treatment to a threshold of 2% and compare the distribution of the effective parameters between the cases with variability set to 0.25 and 0.50, we found that the ranges of values were wider in the second case. Altogether the results presented in this section indicate that the higher the variability in the population, the harder to find effective treatments. Nonetheless, the results also show that even with high values of parameters variability, it is still possible to find very effective treatment strategies.

4 Discussion

In this work, we used a hybrid multi-scale model that merges a mass-action kinetics model of the TNF receptor with a cancer cell Boolean model of different signaling pathways. Moreover, these models are embedded in an agent-based framework that allows considering populations of cells in a defined microenvironment. By performing a model exploration, we have shown that the effective treatments parameter can be found in different cells’ geometries, including 2D monolayers and 3D spheroids. Furthermore, by performing a uniform random sampling of the effective treatment spaces, we found that the parameters for 2D and 3D arrangements exhibit different distributions. These differences are more pronounced in the case of the TNF concentration and the pulse duration, where the effective treatments for the 3D spheroid case are notably more constrained than those found for the 2D disc. We hypothesize that the 3D configuration imposes spatial constraints in the diffusion of the TNF, which restraint the space of values for the candidate’s effective treatment.

We also found that some parameters’ combinations exhibit correlations, indicating that one parameter change can be compensated by adjusting another one. For instance, we observe that the pulse period positively correlates with TNF in the 2D and 3D. This means that increasing the period between pulses can be compensated by increasing the pulse concentration. Interestingly, the correlation between these two parameters is stronger in the 3D spheroid than in the 2D monolayer, showing how the former case is more constrained than the latter. For the 3D spheroid, we also found a strong positive correlation between the pulse period and its duration, showing that these two parameters can also compensate for each other. Nonetheless, the correlation between these two parameters is very low and does not show statistical significance in the 2D monolayer. Finally, we found a negative correlation between the pulse concentration and its duration in the 2D monolayer showing that, in these cases, an increase in the concentration of the pulse can be compensated by a reduction of its duration. Altogether these results indicate that the structure of parameters spaces of effective treatment depends on the spatial cell distribution. Therefore, treatment strategies that work in a 2D monolayer may not work in a 3D spheroid.

Although our model only considers tumor cells, the total drug supplied will be critical when healthy cells are also present. For this reason, we analyzed the total TNF supplied during each experiment as a way to estimate the cytotoxicity associated with the effective treatments. Our results showed that, in general, the distribution of this value tends to be skewed to lower values. We also found that the minimum value of the total TNF for the 2D and 3D are quite similar. Nonetheless, the specific treatment’s parameters are very different; in the 2D case, the effective treatments that use the minimum TNF values has a period around twice times larger than the one in 3D, but the duration of the pulse has half of the duration; the pulse concentration is around four times higher in the 2D case. This analysis also indicates that spatial cell distribution is important for the design of efficient strategies.

We later investigate treatment optimization using two different evolutionary algorithms: GA and CMA-ES. Our results showed that both algorithms could quickly converge to effective treatments, but while the GA can find candidates in the first iteration, the CMA-ES converge to a more robust region of the parameter space. Furthermore, the CMA-ES also finds a statistical distribution for the region of effective treatments. Strikingly, both algorithms converge to slightly different regions of the parameter space in both the 2D and the 3D arrangements. This suggests that the fitness landscape of effective treatments is very rough, exhibiting several valleys of effective treatments regions.

In order to unravel the molecular mechanism behaving like the effective treatments, we analyzed the coarse-grained dynamics of the TNF receptor models for working and none working treatment strategies. The non-working strategies can be grouped into two classes. On the one hand, there are those sets of treatment parameters that do not allow the receptor to reach the threshold needed to propagate the signal downstream of the Boolean model. On the other hand, there are treatment strategies that keep the activation of the receptor for long periods enough to induce cell resistance. In this context, if we define drug resistance as the inability of a cell to respond to a given drug, we find that resistance can come from two aspects: from the dynamics of the Boolean model in response to TNF and from the characteristics of the TNF receptors. The effects of TNF in this Boolean model reported by Calzone et al. (2010) are multifaceted: TNF triggers cells to go from a naive to a survival state but to commit cells to necrosis and apoptosis. Once the cells are committed to either survival, necrosis, or apoptosis, they cannot go back, causing a resistance that can be due to phenotypic variability. As this model was studied using a stochastic Boolean simulator (Stoll et al., 2012), it was possible to capture its dynamics and see that these commitments were not equally fast, or even, that there was a window of activation that allowed controlling the commitment to survival and commiting the cells to necrosis (Letort et al., 2018).

In addition to the Boolean model, our hybrid model also has a mass-action kinetics model that can cause another type of resistance. As we see in Section 3.2 and Section 3.3, there are values for the receptor’s kinetic parameters that prevent the cell from the regulatory effects produced by the binding of TNF. Therefore, when we consider heterogeneous populations by introducing variability in the kinetic parameters of the TNF receptor, we observed that beyond a critical value of the parameter that controls variability, most of the effective treatments that work in the homogeneous population fail to reduce tumor growth. Our hypothesis is that, with high variability, some of the cells could have kinetic parameters that make them insensitive to the treatment, and thus they will produce the relapse after the sensitive ones have been killed by the treatment. We found that this is the main cause of the non-optimal parameters sets found by the optimization techniques within heterogeneous populations.

Regarding the critical value for the control parameter, it is different for the 2D and the 3D cell arrangements, with a lower value in the latter case. We hypothesize that this difference may be due to two factors. The first is because the space of effective treatments strategies is more constrained in the 3D case. The second reason we propose is due to the differences in the total number of cells simulated in the 2D and 3D. While we set the same radius for the disc and the spheroid, the numbers of initial cells are ∼150 and ∼1,000, respectively. Because variability is generated by the sampling of random parameters, a larger number of cells increases the probability of getting a set of kinetic parameters that make the cell insensitive to the TNF. We have shown that population variability can cause resistance. The higher the variability, the harder to find effective treatments. However, even in the cases of the maximum variability analyzed, the algorithms can find a few sets of candidate effective treatments.

5 Conclusion

Multi-scale modeling allows for gaining mechanistic insights in dynamic drug dosages and predicting novel strategies for treatments. Even though, in the last few years, it has been great progress in the field (Montagud et al., 2021), it is known that virtual drug screens seldom match with clinical trials results. Thus, we need to acknowledge that we are far from using these models at the patient’s bedside (Horvath et al., 2016). One of the improvements that would help close this gap would be to have simulations and optimizations that account for and embrace uncertainty. We have hereby presented a free-to-use, open-source framework that allows optimizing treatment strategies with varying levels of uncertainty. We tested the framework using a multi-scale model of cancer growth in different cell arrangements introducing population variability to show that population heterogeneity is critical, either caused by the cells’ state, their parameters, or the population size, affecting the optimal parameter sets. We found that our model exploration workflow can find effective treatments in all the studied conditions. Most importantly, our results show that cells’ spatial geometry and population variability should be considered when optimizing treatment strategies to find robust parameter sets.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found at: https://github.com/bsc-life/spheroid-tnf-v2-emews.

Author Contributions

MP-L and AV conceived the work. MP-L and AM implemented the model. MP-L, CA, TD and JS implemented the workflow and the code for analyzing simulation results. MP-L, CA and AM prepared the manuscript. All authors edited, commented, and agreed on the final version of the manuscript.

Funding

This research has received funding from the Horizon 2020 INFORE Project, GA n° 825070 and the Horizon 2020 PerMedCoE Project, GA n° 951773.

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.

Acknowledgments

The authors acknowledge the technical expertise and assistance provided by the Spanish Supercomputing Network (Red Española de Supercomputación) and the computer resources used: the LaPalma Supercomputer (BCV-2019-2-0008), at the Instituto de Astrofísica de Canarias and MareNostrum4 (BCV-2020-3-0016), located at the Barcelona Supercomputing Center.

Supplementary Material

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

Footnotes

¹https://github.com/bsc-life/spheroid-tnf-v2-emews

References

Akasiadis, C., Ponce-de-Leon, M., Montagud, A., Michelioudakis, E., Atsidakou, A., Alevizos, E., et al. (2021). Parallel Model Exploration for Tumor Treatment Simulations. Computational Intelligence. doi:10.1111/coin.12515

CrossRef Full Text | Google Scholar

An, G. (2010). Closing the Scientific Loop: Bridging Correlation and Causality in the Petaflop Age. Sci. Transl. Med. 2, 41ps34. doi:10.1126/scitranslmed.3000390

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, A. R. A., Weaver, A. M., Cummings, P. T., and Quaranta, V. (2006). Tumor Morphology and Phenotypic Evolution Driven by Selective Pressure from the Microenvironment. Cell 127, 905–915. doi:10.1016/j.cell.2006.09.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Brady, S. W., McQuerry, J. A., Qiao, Y., Piccolo, S. R., Shrestha, G., Jenkins, D. F., et al. (2017). Combating Subclonal Evolution of Resistant Cancer Phenotypes. Nat. Commun. 8, 1231. doi:10.1038/s41467-017-01174-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Calzone, L., Tournier, L., Fourquet, S., Thieffry, D., Zhivotovsky, B., Barillot, E., et al. (2010). Mathematical Modelling of Cell-Fate Decision in Response to Death Receptor Engagement. Plos Comput. Biol. 6, e1000702. doi:10.1371/journal.pcbi.1000702

PubMed Abstract | CrossRef Full Text | Google Scholar

[Dataset] Holland, J. (1975). Adaptation in Natural and Artificial Systems. ann arbor: University of michigan press.

Google Scholar

Ecker, S., Pancaldi, V., Rico, D., and Valencia, A. (2015). Higher Gene Expression Variability in the More Aggressive Subtype of Chronic Lymphocytic Leukemia. Genome Med. 7, 8. doi:10.1186/s13073-014-0125-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Fischer, R., Maier, O., Naumer, M., Krippner-Heidenreich, A., Scheurich, P., and Pfizenmaier, K. (2011). Ligand-induced Internalization of TNF Receptor 2 Mediated by a Di-leucin Motif Is Dispensable for Activation of the NFκB Pathway. Cell Signal. 23, 161–170. doi:10.1016/j.cellsig.2010.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Fitzpatrick, J. M., and Grefenstette, J. J. (1988). Genetic Algorithms in Noisy Environments. Mach Learn. 3, 101–120. doi:10.1007/bf00113893

CrossRef Full Text | Google Scholar

Flobak, Å., Baudot, A., Remy, E., Thommesen, L., Thieffry, D., Kuiper, M., et al. (2015). Discovery of Drug Synergies in Gastric Cancer Cells Predicted by Logical Modeling. Plos Comput. Biol. 11, e1004426. doi:10.1371/journal.pcbi.1004426

PubMed Abstract | CrossRef Full Text | Google Scholar

Fortin, F.-A., De Rainville, F.-M., Gardner, M.-A., Parizeau, M., and Gagné, C. (2012). DEAP: Evolutionary Algorithms Made Easy. J. Machine Learn. Res. 13, 2171–2175. doi:10.5555/2503308.2503311

CrossRef Full Text | Google Scholar

Frieboes, H. B., Edgerton, M. E., Fruehauf, J. P., Rose, F. R. A. J., Worrall, L. K., Gatenby, R. A., et al. (2009). Prediction of Drug Response in Breast Cancer Using Integrative Experimental/Computational Modeling. Cancer Res. 69, 4484–4492. doi:10.1158/0008-5472.CAN-08-3740

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghaffarizadeh, A., Heiland, R., Friedman, S. H., Mumenthaler, S. M., and Macklin, P. (2018). PhysiCell: An Open Source Physics-Based Cell Simulator for 3-D Multicellular Systems. Plos Comput. Biol. 14, e1005991. doi:10.1371/journal.pcbi.1005991

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldman, A., Majumder, B., Dhawan, A., Ravi, S., Goldman, D., Kohandel, M., et al. (2015). Temporally Sequenced Anticancer Drugs Overcome Adaptive Resistance by Targeting a Vulnerable Chemotherapy-Induced Phenotypic Transition. Nat. Commun. 6, 6139. doi:10.1038/ncomms7139

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, N., and Ostermeier, A. (2001). Completely Derandomized Self-Adaptation in Evolution Strategies. Evol. Comput. 9, 159–195. doi:10.1162/106365601750190398

PubMed Abstract | CrossRef Full Text | Google Scholar

Horvath, P., Aulner, N., Bickle, M., Davies, A. M., Nery, E. D., Ebner, D., et al. (2016). Screening Out Irrelevant Cell-Based Models of Disease. Nat. Rev. Drug Discov. 15, 751–769. doi:10.1038/nrd.2016.175

PubMed Abstract | CrossRef Full Text | Google Scholar

Jabs, J., Zickgraf, F. M., Park, J., Wagner, S., Jiang, X., Jechow, K., et al. (2017). Screening Drug Effects in Patient‐derived Cancer Cells Links Organoid Responses to Genome Alterations. Mol. Syst. Biol. 13, 955. doi:10.15252/msb.20177697

PubMed Abstract | CrossRef Full Text | Google Scholar

Kessler, D. A., Austin, R. H., and Levine, H. (2014). Resistance to Chemotherapy: Patient Variability and Cellular Heterogeneity. Cancer Res. 74, 4663–4670. doi:10.1158/0008-5472.CAN-14-0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, E., Kim, J.-Y., Smith, M. A., Haura, E. B., and Anderson, A. R. A. (2018). Cell Signaling Heterogeneity Is Modulated by Both Cell-Intrinsic and -extrinsic Mechanisms: An Integrated Approach to Understanding Targeted Therapy. Plos Biol. 16, e2002930. doi:10.1371/journal.pbio.2002930

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, M. J., Ye, A. S., Gardino, A. K., Heijink, A. M., Sorger, P. K., MacBeath, G., et al. (2012). Sequential Application of Anticancer Drugs Enhances Cell Death by Rewiring Apoptotic Signaling Networks. Cell 149, 780–794. doi:10.1016/j.cell.2012.03.031

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, R. E. C., Qasaimeh, M. A., Xia, X., Juncker, D., and Gaudet, S. (2016). NF-κB Signalling and Cell Fate Decisions in Response to a Short Pulse of Tumour Necrosis Factor. Sci. Rep. 6, 39519. doi:10.1038/srep39519

PubMed Abstract | CrossRef Full Text | Google Scholar

Letort, G., Montagud, A., Stoll, G., Heiland, R., Barillot, E., Macklin, P., et al. (2018). PhysiBoSS: A Multi-Scale Agent-Based Modelling Framework Integrating Physical Dimension and Cell Signalling. Bioinformatics 35, 1188–1196. doi:10.1093/bioinformatics/bty766

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Yin, Q., and Wu, H. (2013). Structural Basis of Signal Transduction in the TNF Receptor Superfamily. Adv. Immunol. 119, 135–153. doi:10.1016/B978-0-12-407707-2.00005-9

PubMed Abstract | CrossRef Full Text | Google Scholar

McGranahan, N., and Swanton, C. (2015). Biological and Therapeutic Impact of Intratumor Heterogeneity in Cancer Evolution. Cancer Cell 27, 15–26. doi:10.1016/j.ccell.2014.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Metzcar, J., Wang, Y., Heiland, R., and Macklin, P. (2019). A Review of Cell-Based Computational Modeling in Cancer Biology. JCO Clin. Cancer Inform. 3, 1–13. doi:10.1200/CCI.18.00069–

PubMed Abstract | CrossRef Full Text | Google Scholar

Montagud, A., Ponce-de-Leon, M., and Valencia, A. (2021). Systems Biology at the Giga-Scale: Large Multiscale Models of Complex, Heterogeneous Multicellular Systems. Curr. Opin. Syst. Biol. 28, 100385. doi:10.1016/j.coisb.2021.100385

CrossRef Full Text | Google Scholar

Ozik, J., Collier, N., Heiland, R., An, G., and Macklin, P. (2019). Learning-accelerated Discovery of Immune-Tumour Interactions. Mol. Syst. Des. Eng. 4, 747–760. doi:10.1039/C9ME00036D

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozik, J., Collier, N. T., Wozniak, J. M., Macal, C. M., and An, G. (2018b). Extreme-scale Dynamic Exploration of a Distributed Agent-Based Model with the Emews Framework. IEEE Trans. Comput. Soc. Syst. 5, 884–895. doi:10.1109/TCSS.2018.2859189

PubMed Abstract | CrossRef Full Text | Google Scholar

Ozik, J., Collier, N. T., Wozniak, J. M., and Spagnuolo, C. (2016). “From Desktop to Large-Scale Model Exploration with Swift/T,” in 2016 Winter Simulation Conference (Washington, DC, United States: WSC), 206–220. doi:10.1109/WSC.2016.7822090

CrossRef Full Text | Google Scholar

Ozik, J., Collier, N., Wozniak, J. M., Macal, C., Cockrell, C., Friedman, S. H., et al. (2018a). High-throughput Cancer Hypothesis Testing with an Integrated PhysiCell-EMEWS Workflow. BMC Bioinformatics 19, 483. doi:10.1186/s12859-018-2510-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Reuillon, R., Leclaire, M., and Rey-Coyrehourcq, S. (2013). Openmole, a Workflow Engine Specifically Tailored for the Distributed Exploration of Simulation Models. Future Generation Comput. Syst. 29, 1981–1990. doi:10.1016/j.future.2013.05.003

CrossRef Full Text | Google Scholar

Sedger, L. M., and McDermott, M. F. (2014). TNF and TNF-Receptors: From Mediators of Cell Death and Inflammation to Therapeutic Giants - Past, Present and Future. Cytokine Growth Factor. Rev. 25, 453–472. doi:10.1016/j.cytogfr.2014.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaffer, S. M., Dunagin, M. C., Torborg, S. R., Torre, E. A., Emert, B., Krepler, C., et al. (2017). Rare Cell Variability and Drug-Induced Reprogramming as a Mode of Cancer Drug Resistance. Nature 546, 431–435. doi:10.1038/nature22794

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoll, G., Viara, E., Barillot, E., and Calzone, L. (2012). Continuous Time Boolean Modeling for Biological Signaling: Application of gillespie Algorithm. BMC Syst. Biol. 6, 116. doi:10.1186/1752-0509-6-116

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, X., Bao, J., and Shao, Y. (2016). Mathematical Modeling of Therapy-Induced Cancer Drug Resistance: Connecting Cancer Mechanisms to Population Survival Rates. Sci. Rep. 6, 22498. doi:10.1038/srep22498

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, K. S., Man, K. F., Kwong, S., and He, Q. (1996). Genetic Algorithms and Their Applications. IEEE Signal. Process. Mag. 13, 22–37. doi:10.1109/79.543973

CrossRef Full Text | Google Scholar

Whitley, D. (1994). A Genetic Algorithm Tutorial. Stat. Comput. 4, 65–85. doi:10.1007/bf00175354

CrossRef Full Text | Google Scholar