Modeling small cell lung cancer (SCLC) biology through deterministic and stochastic mathematical models

Abstract

Mathematical cancer models are immensely powerful tools that are based in part on the fractal nature of biological structures, such as the geometry of the lung. Cancers of the lung provide an opportune model to develop and apply algorithms that capture changes and disease phenotypes. We reviewed mathematical models that have been developed for biological sciences and applied them in the context of small cell lung cancer (SCLC) growth, mutational heterogeneity, and mechanisms of metastasis. The ultimate goal is to develop the stochastic and deterministic nature of this disease, to link this comprehensive set of tools back to its fractalness and to provide a platform for accurate biomarker development. These techniques may be particularly useful in the context of drug development research, such as combination with existing omics approaches. The integration of these tools will be important to further understand the biology of SCLC and ultimately develop novel therapeutics.

KEY POINTS

Novel therapies are urgently needed for patients with small cell lung cancer (SCLC)

The use of mathematical models in medicine is becoming progressively more robust due to the exponentially increasing processing power of new hardware and software

Computational models thus far have mostly focused on comprehensive cancer biology models rather than individualized disease models

Mathematical modeling integrated with computational modeling can be used to simulate tumor growth, mutational heterogeneity, cellular automata, and mechanisms of metastasis

Mathematical oncogenic models can provide insight into biological mechanisms of action, identification of biological markers, quantification of image analysis and understanding of chaotic cell dynamics, especially in immune function

MAIN

Self-similar, repeating patterns, also known as fractals, can describe the universe and mimic nature’s highly complex structures. Nature often presents itself in fragments of patterns that are similar, yet not identical. The conservation of structural patterns in nature is manifested across organisms from different realms. One prominent example is the respiratory analogy between trees and lungs, where conservation of both fractal design (self-similarity) and function is elegantly manifested (Figure 1). The central tree trunk/trachea divides into wider branches/bronchi, which bifurcate into increasingly smaller branches/bronchioles, and eventually conclude in leaves/alveoli. The synchrony in morphology is also translated into function, where gas exchange occurs at the same terminal end, leaves or alveoli. Maintenance of this stable, yet complex natural order controls biological equilibrium and thus life. Quantifying changes in chaotic patterns could provide a valuable diagnostic tool to distinguish at an early stage between healthy and abnormal tissue, paving the way for more effective and personalized cancer treatments.

Figure 1: Fractal geometry in SCLC. (A) 3D visualization of SCLC patient lungs and their fractal properties. (B) Fractal tree generated using Swift 3.1 programming language. (C) Histology of SCLC showing the fractal geometry of tumor tissue. Hematoxylin and eosin stained SCLC specimen. SCLC tends to travel in clusters and the encircled areas are examples of this clustered tumor property. Specifically, a cluster infiltrating into a blood vessel is shown on the left.

Our goal is to develop an understanding of the mathematical basis of cancer in the context of small cell lung cancer (SCLC) growth, mutational heterogeneity, and mechanisms of metastasis. The results can serve as a biomimetic framework for organizing and applying the wealth of available genomic and clinical data. We aim to communicate this mathematical understanding to clinicians and tumor biologists in the hope of generating precise treatment maps and therapies with improved outcomes. It is hypothesized that cancer is switched on by a chaotic imbalance through stochastic processes of mutation generation and genetic drift, which are inherently random. However, further malignant development seems to recover this balance, albeit to one far from normal, in favor of a deterministic process through clonal selection. Deterministic models include ordinary and partial differential equations that describe tumor kinetics and stroma, whereas stochastic models include cancer initiation and progression such as cellular automata and group theory. These tools may provide a focal point for integrated systems approaches that are aimed at the development of new cancer therapeutics.

SMALL CELL LUNG CANCER (SCLC) AND ITS CLINICAL STATUS

SCLC is thought to originate from neuroendocrine cells that normally reside in the lung epithelium [1], still retaining some of the neuroendocrine cell markers, including CD56 [2, 3], chromagranin, and synaptophysin [1] or the neuroendocrine transcription factor ASCL1 (achaete-scute complex homolog-like 1 (ASCL1)) [4]. SCLC is typically highly metastatic in the majority of patients with a large number of circulating tumor cells (CTCs) in the periphery, likely contributing to relapse and poor prognosis [5, 6]. CTCs can have stem cell properties and may help detect early disease, define treatment efficacy or serve as diagnostic markers using liquid biopsy methodology [7–9]. The majority of SCLC patients have a history of smoking and the disease can present itself many years after smoking cessation [10, 11]. SCLC is characterized by rapid and early metastatic growth. SCLC itself is typically not confined to a single tumor mass, instead it spreads rapidly, often as clusters of cancer cells. These clusters or spheroids are commonly observed during metastasis and may contribute to chemoresistance [12]. This suggests that molecules that regulate cell/cell interactions, cell/matrix interactions or in general cytoskeletal functions, are potential therapeutic targets for SCLC.

SCLC is either staged as limited disease (LD), with tumors localized to one hemithorax with potential regional lymph node involvement or ipsilateral pleural effusion, or it is staged as extensive disease (ED) with expanded metastatic states [13]. Most patients are not diagnosed at an early stage, have metastatic disease and face limited treatment options with an overall 5-year survival of less than 7% [14]. It seems apparent that early diagnosis using appropriate biomarkers would improve outcomes. However, to date there are no effective methods for early screening and detection for SCLC. With a 5-year survival rate of less than 20% and more than 30,000 deaths per year, SCLC is defined as a “recalcitrant” cancer in the US. Treatment of SCLC varies depending on the extent of the disease. Standard therapy for LD includes chemotherapy, mostly in combination with radiation. Currently, approximately a quarter of patients with LD can be cured with current standard treatments. Typically, platinum-based chemotherapeutics are combined with etoposide and once or twice daily thoracic radiation therapy before the third cycle of chemotherapy [15–18]. In general, SCLC initially responds well to chemotherapy, yet eventually resistance develops. The median survival for patients diagnosed with LD is 15–20 months and 5.5 months for patients diagnosed with ED. Chemotherapy is also the first-line treatment for patients with ED, but treatment is often palliative [19]. For patients that are not sensitive to etoposide/carboplatin, second-line treatment is mainly limited to the topoisomerase I inhibitor topotecan [20, 21]. Additional therapeutic approaches and targeted therapies have been developed and tested or are under development but have not yet led to a breakthrough in treatment and 5-year overall survival has only marginally improved in the past 40 years [14, 22]. There is thus a clear need to discover new targeted therapies for SCLC and it will be crucial to define biological mechanisms that cause cancer promotion, progression or metastasis and link them to genetic changes. Advances in understanding the biology of SCLC have largely relied on in vitro studies using cell line models and may not be suitable for unraveling the sequence of events leading to the aberrations responsible for tumor initiation. Solving these problems remains a challenge, and clinical application of acquired discoveries is often lacking [23, 24].

SCLC CELLS REPRESENT A DIVERSE POPULATION

Intrinsic and persistent genomic instability in SCLC are a driving force for clonal diversity and disease evolution. With about 175 mutations per tumor, the rate of genomic alterations in SCLC is amongst the highest in solid tumors (5.5 to 7.4 mutations per Mb) [25, 26]. Frequent G-to-T transversions indicate a tobacco carcinogenesis signature and are consistent with a smoking history [25]. Even though SCLC is thought to be of neuroendocrine origin, the cancer stem cell (CSC) population is not well defined and appears diverse on a phenotypical and genomic level, which may be due to high rate of mutations and/or epigenetic regulation. Phenotypically, SCLC CSCs are not well defined and may contain several markers including SOX2, CD44, CD56 (NCAM), CD90, CD105, CD133, Sall4, Oct4, nestin, S100β or vimentin [27, 28]. Interestingly, as the ASCL1 transcription factor regulates neuroendocrine features, it also cooperates with Notch signaling in normal airway stem cell differentiation [29]. In SCLC, this pathway is altered and overexpression of E2F3 as well as loss of RB function drive disease progression, likely supported by loss of function mutations in the tumor suppressor TP53 [30]. Indeed, TP53 (100%) and RB1 (93%) are the highest mutated genes in SCLC cases without chromotripsis. [31]. Further, expression of the histone-lysine methyltransferase enhancer of zeste homolog 2 (EZH2) strongly correlated with disruption of E2F transcription factors/RB1 pathway, found in 96% of SCLC [32]. Changes in target gene profiles are also associated with a “stem cell-like” phenotype [32] and in some cases with acquired chemoresistance [33]. NOTCH family genes are frequent targets of inactivating mutations in SCLC, with about 25% of tumors being affected and this inactivation is required for optimal tumor growth [31]. In addition, expression of Delta-like protein 3 (DLL3), an inhibitory ligand of NOTCH that is regulated through ASCL1, is expressed in more than 80% of SCLC. Inhibition of DLL3 with rovalpituzumab, a DLL3-targeted antibody-drug conjugate, in recurrent SCLC shows single-agent anti-tumor activity and would be expected to also specifically target CSCs [34]. The overall NOTCH pathway mutational profile in SCLC is consistent with suppression of Notch activity [31]. Additional genes that are mutated at a significant level include the G protein-coupled receptor FPR1, the G protein regulatory protein RGS7 as well as KIAA1211 and COL22A1 [31]. There is also amplification of MYC family members in about 16% of SCLC [25], which function through transcriptional regulation [35]. There are also alterations in receptor tyrosine kinases, including FGFR1 and KIT or less frequent targets, including MET, RON and EPHB4 as well as proximal downstream effectors, such as PIK3CA or PTEN [31, 36–38]. Also, these receptor tyrosine kinases may signal in part through reactive oxygen species (ROS) as a result of metabolic reprogramming that can non-specifically activate signaling cascades and contribute to DNA damage [39]. Metabolic changes in cancer cells are intrinsically associated with transformation and may not only provide energy, but also intermediates for anabolic pathways or metabolites that affect gene expression and ultimately result in changes within the tumor microenvironment, leading to an overall growth advantage. There are additional mechanisms that reflect the genetic heterogeneity within this cancer [25, 31].

SCLC METASTASIS AND THE TUMOR MICROENVIRONMENT

The process of metastasis is divided into several characteristic phases, starting with invasion of tumor cells into surrounding tissues and then eventually the blood stream. As circulating tumor cells (CTCs) they can reach distant sites and grow if they have acquired the capability to survive and interact with the various tissues, such as extravasation through the endothelial lining of blood vessels [40]. The success of this process depends on specific cellular properties that may vary and are not only determined by genomic changes and their cellular consequences but also by the effects that the tumor has on its microenvironment and how it interacts or responds to them. Both extravasion and subsequent intravasion to establish metastatic sites, can be regulated at multiple levels, involving ligands within the extracellular matrix ECM, their receptors, including selectins, integrins, cadherins, CD44 and others, or chemokines and cytokines and their receptors. Additional interaction with immune cells or stromal cells further determine metastatic function [41]. A retrospective study analyzed 251 SCLC patients diagnosed between 1999 and 2000 and found 152 (60.6%) with distant metastases. Target organ involvement included 20.3% liver, 18.3% bone, 15.5% brain, 10.0% lung and 6.0% of adrenal gland [42]. (Figure 2; Supplementary Video 1) The model was generated by javascript bubble chart, where the initial state is fixed with all cancer cells positioned at the primary site, and the final state is quantitatively fixed based on the metastasis sites population reported in Nakazawa et al. [42]. Cells with different metastatic phenotypes first appear at the primary site, then cells with similar metastatic phenotypes cluster together and eventually metastasize to the different sites. Metastasis is a multi-step process in the cancer model, we think about mitogenesis, morphogenesis, and motogenesis [43]. The mitogenesis gives the proliferation that also then effects morpho- and moto-genesis. The process of metastasis is dependent on genetic regulation, protein network, and tumor-stroma interaction. As an example, we have shown that PAX5 transcription factor is highly expressed in SCLC [44]. This then regulates chemokine receptor CXCR4 and RTKs such as MET and RON. The receptors themselves cause a plethora of signal transduction events such as activation of the focal adhesion protein FAK and Actin cytoskeleton. Ultimately this leads to increased motility, invasion, and metastasis.

Figure 2: Most common metastasis sites of SCLC. Simulation of early time points and the endpoint (from top left to bottom right) of tumor cells (purple) metastasizing to different organs (colored as indicated). See also Supplementary Video 1.

EMT is a crucial process for metastasis formation in cancer progression and development. A few studies have identified small cell lung cancer as an epithelial to mesenchymal transition (EMT)-like cancer [45, 46]. Small cell lung cancer and other carcinomas that undergo EMT have been observed to have an enhanced metastatic capacity, but this is complicated by the presence of epithelial differentiation at the metastatic sites [47]. The forward and backward transitioning between epithelial and mesenchymal phenotypes has become more well understood and it is believed that cancer cells with hybrid phenotypes, in which the EMT network is in a hybrid oscillating state, are associated with more aggressive cancer behaviors [48–50]. In the role of metastasis, it is theorized that the presence of a high number of hybrid E/M cells may form CTC clusters, which form much more metastases than single CTCs [51–53].

The significance of the tumor stromal microenvironment for tumorigenicity and metastasis in SCLC has been well established and some of the mechanisms involved have been already unraveled. In particular, high levels of extracellular matrix (ECM) proteins are found to be specifically associated with SCLC tumors. ECM proteins bind to integrin and other cell surface receptors, effectively regulating cellular transformation processes through ‘outside-in’ signaling. Conversely, ‘inside-out’ signaling regulated by tumor cells defines the interaction of receptors through changes in their functional interaction with ECM proteins [54, 55]. SCLC can actively participate in producing ECM proteins, such as fibronectin and laminin to support anti-apoptotic mechanisms and reduce the efficacy of chemotherapy [56]. Binding of fibronectin to the β1 integrin receptor stimulate PI3K signaling mechanisms that regulate apoptosis, cell motility/morphology and adhesion [56, 57]. Laminin also signals through the PI3K pathway and this activity may be required for suppression of apoptosis, morphological changes and chemoresistance [58]. Invasion and metastasis requires protease activity to penetrate the barrier provided by the ECM. Expression of the ADAM-12 (A disintegrin and metalloprotease-12) protease was found in 73% of tumor specimen and other ADAM family members in 10-40% of specimen. Targeted knockdown of ADAM-12 in the H1688 SCLC cell line reduced cell growth, invasion and metastatic function [59]. Hyaluronic acid, another ECM component, binds mostly to the CD44 receptor and can be cleaved by hyaluronidase. In SCLC, the expression of CD44 and hyaluronidase was found to be low or absent, thus altering potential metastatic function during intravasion and extravasion [60, 61].

The MET receptor, and its ligand, HGF, may also play an important role in metastatic process as well as allowing cancer cells to survive in distant sites [62]. The juxtamembrane R988C and T1010I mutations in MET introduced in a pre-clinical cell line model can lead to growth factor-independent cell growth. Moreover, when overexpressed in the H446 SCLC cell line, both mutations were sufficient to alter cell morphology, adhesion, foci-formation and soft-agar colony formation [38]. MET activation leads in part to phosphorylation and activation of FAK (focal adhesion kinase), AKT, the ERK1/2 pathway and others. Expression of active MET can be targeted by the MET inhibitor, SU11274 or RNA interference, thus regulating MET-dependent invasion [63, 64]. In vitro testing of metastatic cells in an orthotropic transplant model showed HGF/MET-dependent motility and invasion of metastatic cells and metastasis itself in bone marrow, kidney and brain could be suppressed by MET inhibition [65]. Ligation of the related receptor tyrosine kinase RON by its ligand macrophage stimulating protein (MSP) suggests a rather specific role in promoting liver metastasis in a SCLC cell line-based in vivo mouse model [66]. In patients, there is high expression of the SDF-1/CXCL12 chemokine receptor CXCR4. Normally, SDF-1 is produced in the bone marrow by stromal cells and acts as a homing factor for hematopoietic cells. In SCLC, SDF1 enhances the binding of VCAM-1 (vascular cell adhesion molecule 1 or CD106) to VLA-4 (Very Late Antigen-4 or integrin α4β1). This may explain in part the high propensity of SCLC to metastasize to the bone marrow [67]. This is consistent with an orthotopic xenograph in vivo mouse models, wherein tumor growth could be reduced by cisplatin and etoposide, but formation of metastases was not affected. However, the CXCR4 inhibitor AMD3100, was somewhat less effective in reducing tumor growth but significantly reduced suppressed metastasis formation by 43% [68].

The phenotypic plasticity of small cell lung cancer is a dynamic chaotic system where a variety of starting conditions exhibited by molecular interactions and environmental fluctuations evolves into a set of values of the variables, which can be understood as an “attractor” (steady state) where slight perturbations can have robust outcomes [69, 70]. At the core of this system is an endogenous molecular-cellular network that is regulated by gene regulatory networks (GRNs) where fluctuations from normal attractor to cancer attractor are responsible for cancer formation and the eventual cancer progression as cancer attractors transition phenotypes responsible for resistance [70, 71]. A recent study explored the stability of the cancer state that relies on the cancer attractor and the results implied that at a certain state the cancer can no longer be reversed to a non-malignant phenotype [72]. However, existing data shows that it is possible to reverse and maintain aggressive “cancer attractors” by applying a stepwise therapeutic approach by targeting the dynamic system responsible for the primary factors of metastasis and progression [70]. In one example, the transition from normal state occurs through the stepwise process of MDM2 on, CDK2 on, RB off leading to cancer state [72]. However, the reversal processes is similar where RB goes on, CDK2 and CDK4 off, and off of MDM2 leading to a normal state in that sequence [72]. This method of reversing the cancer attractor is vital for SCLC due to the vital role of RB in SCLC oncogenesis, and the restoration of this tumor repressor gene may be a potential therapeutic tactic [72].

As can be appreciated, the genetics/proteomics and growth characteristics for SCLC are unique. In spite of decades of research, we have not made progress in the therapy for SCLC. The role of the immune system in the biology of SCLC is just beginning to emerge and typically depends on the functional interaction of lymphocytes, tumor cells and antigen presenting cells. Currently, several trials determine the efficacy of immune checkpoint inhibitors in SCLC. This treatment is done with a curative intent, since it has the potential to target cancer stem cells, rather than blocking metastatic spread of the tumor [73]. Some of the additional advances for this disease will come from early detection, understanding biology, developing novel therapies, and prevention of recurrence. In order to understand the various genetic and cellular biology of SCLC, we have started to develop mathematical models. In particular, one can evaluate the growth characteristics with first order kinetics such as ordinary differential equations, tumor-stroma interactions with partial differential equations, and geometry/cell movement/metastasis with chaos theory/fractal analysis.

COMPUTATIONAL MODELING

Discrete models – cellular automata and agent based models

Single cell events such as mutations have been modeled by discrete methods such as cellular automata [74], agent-based models and hybrid continuum-discrete approaches [75]. These discrete models led to multiscale modeling of cancer where by ordinary differential equations (ODEs) are linked to cellular level parameters [76].

The cellular Potts model [77] is a more generalized cellular automata (CA) that uses lattice dynamics to study interactions among biological cells. The cellular Potts model was used to study the formation of cell clusters as a consequence of assuming configurations of minimal adhesive free energy [78]. The CA model and the cellular Potts model fall under agent based models (ABM). The variables in ABM are individuals. These individuals are considered as agents and a set of prescribed rules govern the behavior of these agents. In the case of CA, the individual lattice cells are the agents.

A discrete agent based spatio-temporal model that incorporates the effects of nutrient supply, mechanical confinement that represents the tissue resistance against tumor cell movement and toxicity of metabolites in the context of brain tumor progression was developed by Mansury et al. [79]. They simulated the complex dynamic self-organizing and adaptive processes observed in tumors, namely spatial aggregation of tumor cells as clusters and their migration in search of suitable survival conditions.

Hybrid agent based models combine the unusual effectiveness of continuum deterministic models to capture tumor dynamics at the tissue scale with discrete CA models at cellular and subcellular scales [80]. Tumor invasion of stroma and surrounding tissue are modelled as coupled non-linear partial differential equations (PDEs). The PDEs are discretized to model cell migration and form the basis of the hybrid discrete-continuum model. This model enables specific properties of cells to be described such as proliferation, death, cell-cell adhesion and mutation.

The software BioFVM [81] was utilized to simulate the cellular automaton models with some adjustments to values that would favor SCLC (Figure 3, Supplementary Videos 2 and 3). To do this, by increasing the birth rate by .025 and increasing the viable cell Hill coefficient to 4 to simulate our growth rate in Figure 4, modified the viable life span of the tumor cells to increase the probability of apoptosis over time and decrease the time between apoptosis and necrosis due to hypoxia, and introduced a 3% rate of necrosis for perinecrotic tumor cells to fit our carrying capacity obtained in Figure 4. The remaining parameters of the simulation remain similar, where the duration of necrosis and cell death are the same, the system begins with a single cell at time = 0, and the drug parameters are unchanged. The first action performed by the BioFVM program as it creates the simulation is to fill any unoccupied spaces with fluid and considers any dead cells as an empty space [81]. The model then uses the following reaction-diffusion equations for the oxygen and drug, respectively, to apply the uptake rate across each cell:

Figure 3: The cellular automaton model of SCLC growth and necrosis. (A) Growth of the tumor (scale on left side) and oxygen levels within the tumor (color scale on right side) at different time points. See also Supplementary Video 2. (B) Growth of the tumor showing the extent of the necrosis at different times (red = live cells, green = apoptotic cells, and blue = necrotic cells). The necrosis of the tumor initiates at > 180 hours. See also Supplementary Video 3.

Figure 4: Tumor growth models. (A) Tumor growth with dynamic carrying capacity and metastatic burden. The growth of the tumor’s cell population c (blue) and carrying capacity K (black) with λ = 0.192, φ = 5.85, and ϕ = 0.00873, parameters used by Enderling and Chaplain [90]. The evolution of metastatic growth rate d (red) is presented with the parameters from Benzekry et al. [114]. (B) Histogram of tumor evolution (green) alongside the histogram (red) of the metastatic growth. As the primary tumor reaches the maximum carrying capacity, the metastatic burden will increase exponentially from the detached primary cells traveling and growing at distant sites (see curve). The metastatic cell population grows beyond the original carrying capacity until it reaches a new carrying capacity due to the model describing the growth of cells in all metastatic sites rather than just the primary site.

where the treatment by the drug is set at 5 μM at time t = 528 (Day 22). When observing the two-dimensional image of the tumor over time, the cells become more hypoxic the closer they are to the middle. The oxygen concentration in the tumor begins to decrease in steps as the tumor begins to grow, leaving the population of hypoxic cells to grow rapidly.

With the presence of the drug at time t = 528, the live tumor cell’s exposure to the drug E and its response to the drug R (having a Hill coefficient h = 1) are given by the equations below:

where α is the exposure for a half-maximum effect [81].

The BioFVM simulations may be captured after the 2-Dimensional tumor profile, cell automata model, basic agent model, and MATLAB model scripts have been edited and saved and the simulations have been run through the command prompt [75]. Once the simulations have finished running, MATLAB is automatically prompted and a series of images, or visual interpretations of the system at each recorded time interval, are then opened. These images may then be saved as image files. This spatio-temporal model could be utilized to predict and personalize patient response to drug therapy using organoids, spheroids, mouse models, and zebrafish models.

Discrete models – cellular automata, tumor microenvironment, and cell viability

Cancer models based on differential equations address a continuum of cells at the tissue scale where the effect of individual cells is averaged. On the other hand, discrete models of tumor growth based on CA capture the response of individual cells as they interact with one another as well as with its microenvironment. CA is a collection of cells arranged on a lattice that evolve with time according to a defined set of rules that includes the values of neighboring cells. The first practical application of CA was shown by John Conway in “Game of Life” [82]. CA models have been used to study tumor growth at the cellular scale [83–87] and subcellular scale [88, 89]. For example, CA was used to study the invasion of cancerous cells in a population of normal cells by Qi et al. [86]. In this context, a lattice cell represented a single biological cell. The state of each of the cells of the CA was assumed to be normal, cancerous, complex (cancerous and bound by white blood cell) or dead cancerous. Probabilistic rules were then applied to study the dynamics of the cellular states. However, this model did not explicitly consider growth promoting factors (such as presence of blood vessels, nutrient supply and oxygen) and growth inhibiting factors (such as toxic metabolites) for tumors that ‘motivate’ them to move far away from primary sites.

The microenvironment-dependent birth rate b_i and death rate d_i are then modeled how Macklin et al. described as shown below:

where d^* _i,N is a constant [81]. The microenvironment apoptosis rate d_i,A is modeled as:

where d^max is the maximum rate of apoptosis [81].

For the time interval [t, t+Δt], each viable tumor cell has the chance to divide, undergo apoptotic death, or reach necrotic death. The probability for a live cell to perform one of the three actions in each time step until death is described by the equations below [81]:

Each dead cell then has the probability that it will become lytic and rupture, leaving an empty space in the automaton model, based on the average duration the cell death. The probability of the dead cell to reach lysis is expressed in the equation below where T_D is the duration of cell death [81]:

Continuous models − ODEs and PDEs

For a given set of initial conditions, models that produce the same results each time they are solved are known as deterministic models. These differ from stochastic or probabilistic models in that the model results change each time they are solved even though the initial conditions don’t. Deterministic models with one independent variable (“time”) and one or more dependent variables (such as “substrates” or “metabolites”) and represented by ODEs are ideal to capture dynamical processes. For example, Enderling and Chaplain [90] studied the rate of tumor growth cells. Utilizing parameters α (fraction of dividing tumor cells) and β (fraction of tumor dying cells), they showed that tumor cells could either be in 1) quiescence (α – β = 0), 2) proliferative (α > β) or 3) depleting (α < β). Since tumors do not grow indefinitely in size, a more realistic representation of rate of tumor growth should take into account the carrying capacity constraint, K, representing the maximum population of cells also known as carrying capacity of the host cell. Hahnefeldt and colleagues [91] modelled K as a function of time and tumor size as follows:

where ϕ and φ represent constant positive rates of angiogenesis stimulation and inhibition, respectively (Figures 4 and 5, Supplementary Video 4).

Figure 5: Simulation of tumor heterogeneity and tumor burden in SCLC. (A) Tumor heterogeneity and cellular automata at different time points. (B) Tumor burden at corresponding time points, indicating the growth of wild-type cells (red) and cells containing mutations (blue). See also Supplementary Video 4.

Mathematical models solely based on ODEs describe the total number of tumor cells over time but do not consider any spatial variables. It is essential to model the spatial variables along with time since processes such as cancer invasion and metastases are more potent killers than local tumor growth and are inherently spatial in nature. Models based on PDEs such as reaction diffusion are apt for quantitative substances of interest in cancer modeling (such as nutrients or oxygen) at a specific position (space) and time (t). PDE based models are also referred to as continuum models since they are solved for continuously in space and time variables.

For example, Gatenby and Gawlinski [92] were one of the earliest to model cancer invasion as a spatio-temporal evolution of tumor cells (C), enzymes with H+ ions (m) and extracellular matrix (υ) as follows:

where D_C is the diffusion coefficient constant, ρ is the tumor cell proliferation rate constant, δ is the H+ ions production and decay rate constant and γ is the extracellular matrix degradation rate. Appropriate initial conditions and spatial region need to be specified to solve where tumor cells are assumed to proliferate and undergo nonlinear diffusion and secrete H+ ions which diffuse and degrade the normal tissue. The H+ ions are assumed to undergo linear decay with logistic growth for normal tissue in the absence of any cancer cells. Cancer cell migration processes such as haptotaxis (i.e. directional cell migration in response to gradients of cellular adhesion molecules in the extracellular matrix or gradients of the extracellular matrix density) were modeled using a modified version of by Anderson et al. [93].

The PDE based models can be discretized using finite-difference approximations. To study individual cell movement, Anderson et al. investigated the discrete form of the continuous version built to study haptotaxis [94]. Spatial variables were discretized retaining time, t, to be continuous. Stochastic movement rules were incorporated to derive a biased random walk governing the motion of a single tumor cell. Dynamical models of cancer growth leading to chaotic behavior have also been reported [95]. Itik and Banks were able to explicitly show the existence of deterministic chaotic dynamics by modeling the interactions and competitions between tumor cells and other cells of the body such as healthy host cells and activated immune system cells. Based on ideas from Lie algebra [96], the control of chaotic dynamics of cancer growth has been recently formulated [97] in a three dimension cancer model (TDCM) for tumor growth. This spatio-temporal heterogeneity model could be utilized to understand the tumor evolution over time as well as attempt to predict the genetic phenotype that may correlate with metastasis and cancer progression. As we go forward from here, we have to be able to incorporate mathematical modeling in SCLC. We can envision the utilization of the various models in the behavior of cell lines/three-dimensional models, organoids/spheroids along with PDX/CDX models, as well as tumor behavior in natural progression and/or therapeutic response. In particular, we should be able to study the potential for mechanisms of resistance.

MultiFractals – applications to image data analysis

Most naturally occurring phenomena exhibit complexity and irregularity in structure. Classical techniques in mathematics are not adept at capturing the underlying geometry. Irregular objects for which the “Hausdorff dimension” is strictly greater than the “Topological dimension” were defined as fractals by Benoit B. Mandelbrot [98]. Even though fractal objects seem to have complex shapes and exhibit dynamic behavior, they always exhibit a nesting of statistically self-similar features at different scales [99]. The branch of a tree has similar structure to the tree but at a smaller scale. A high degree of self-similarity has also been documented in human tissue [100]. Fractals are compared by their fractal dimension, a real number, which is a measure of the size of the underlying fractal sets. There are multiple methods to compute the fractal dimension. The similarity dimension [101], D_S is ideal for computing the fractal dimension of artificially generated fractals (following precise creation rules and having regular structure). It is defined as follows:

where N represents the number of non-overlapping copies of the object and r (< 1) represents the scaling ratio. The box counting dimension [102] that is most widely used in practice to compute the fractal dimension takes into account the relationship between non-empty boxes and box-size necessary to cover the fractal object. The fractal dimensions of most objects found in nature are best computed using the box counting method. The fractal dimension has been applied as one of the features for classifying pathological tissue in mammograms and tumor blood vessels [103, 104].

The box counting approach to computing the fractal dimension assumes that the fractal object can completely be described by two states (presence inside or outside the box). However, most natural phenomena such as greyscale image of tissue or cells require the intensity inside a box to be considered as well. Since the intensity levels of such an image can be non-uniform resulting in low to high densities, their structure characterized by their fractal dimension would vary with the observed scale.

Unlike synthetically generated fractals, the structure of naturally occurring fractal objects at different scales will be similar but not exactly same as the whole. This results in the co-existence of multiple fractal subsets with different scaling behavior. Therefore, instead of one set of fractals we are typically faced with multifractals.

Multifractal analysis requires the description of a local and global measure over a region. For each point of the region, the local singularity coefficient (i.e. the local dimension) is known as the “Holder exponent” or α values [101, 102]. In image data analysis, α values reflect the local behavior of the intensity measure that are of 4 different types [105, 106]. The corresponding image is often referred to as the α –image. A region would contain a range of positive, finite α values, with a minimum value, α_min and a maximum value, α_max. The fractal dimension can be computed using the box counting method [107] over a set of points (isolated or otherwise) having the same α value. A plot of fractal dimension against α values results in a multifractal spectrum, f(α) (i.e. the global measure) after adjusting for noise in digital images.

Multifractal theory has been applied in the context of medical image data analysis [103] and object classification [101]. Using multifractal analysis of breast cancer tissue prior to chemotherapy Vasiljevic et al. [107] showed that the tissue can be differentiated based on their sensitivity to chemotherapy. Small cell lung cancer embodies the essence of the aberrant tumor cell which despite its diverse genotype has a distinctive, discrete, and easily recognizable tumor tissue [69]. These unique spatial characteristics could be characterized utilizing the multifractal model at different scaling magnifications such as DNA random walks, tumor tissue, and radiological scans.

CONCLUSIONS

We have summarized and applied several mathematical models to SCLC that enable us to replicate some of its unique fractal characteristics. Features of SCLC tumors are a reflection of alterations in genes associated with the disease including overexpression, somatic mutation and amplification, resulting in a plethora of targeted therapies existing or in development including FAK inhibitors, RTK inhibitors targeting KIT, IGF-1R, EPHB4, RON and MET as well as their downstream pathways including the PI3K/AKT/mTOR and MYC [108]. Further areas for drug therapies in SCLC include the apoptotic pathway, the hedgehog and DNA repair pathways, heat shock proteins and HDAC as well as angiogenesis pathways and others [109]. Additional novel therapies in development include targeting of cancer stem cells and immunotherapy [110]. Despite these targeted therapies, the prognosis for SCLC remains grim [111, 112]. Most recently, there have been immunotherapies that have come to fruition. Checkpoint inhibition appears to work with pembrolizumab or nivolumab/ipilimumab. There are specific other programs, such as bispecific antibodies and CAR-T cells that are just beginning to be employed. The effective utilization of mathematical models to these therapies is hampered in part by the availability of quantitative in vitro experimental data and clinical results. Whereas the goal remains to link outcomes to the fractalness of SCLC, these techniques may be particularly useful in the context of drug development research, particularly in combination with existing research platforms, such as various omics approaches that provide large amounts of data. As clinical results can vary widely, dependent on the patient’s specific tumor alterations, it will be important to link them to these genetic changes. Even though the models presented here are not perfect, they provide us with the tools to allow us to breakdown the mechanisms of growth and metastasis and eventually add additional parameters. We can really state that models depend on the circumstance. ODE can be utilized to study tumor kinetics, PDE on tumor-stroma interactions, and chaos theory on geometry and symmetry breaking [113] in metastasis. The complexity of any physical or biological system can be quantified in terms of combinatorial, geometric and functional components. Each of these components is best characterized by their symmetries. Symmetric breaking occurs when the natural symmetry is broken either explicitly or spontaneously. Understanding cancer processes like metastasis as a sequence of symmetry breaking events will be required. The fractal analysis gives us a unique insight on geometry and ultimately could serve as a biomarker. An important goal is to link the existing comprehensive set of analytic tools back to the fractalness of cancer and provide a platform for accurate biomarker development. The identification of appropriate biomarkers in combination with mathematical modeling has the potential to provide breakthroughs in the development of therapeutics. Since SCLC is a rapidly growing disease, it will be important to utilize these powerful mathematical tools to enhance our understanding of the biology of this devastating cancer and ultimately accelerating the development of novel therapeutics.

Abbreviations

SCLC: small cell lung cancer; ASCL1: achaete-scute complex homolog-like 1; LD: limited disease; ED: extensive disease; Mb: Megabase; CSC: cancer stem cell; SOX2: Sex Determining Region Y-Box 2; Sall4: Spalt like transcription factor 4; Oct4: Octamer-binding transcription factor 4; E2F3: E2F transcription factor 3; TP53: Tumor protein p53; RB1: RB transcriptional compressor 1; NCAM: Neural cell adhesion molecule or CD56; EZH2: enhancer of zeste homolog 2; DLL3: Delta-like protein 3; FPR1: formyl peptide receptor 1; RGS7: regulator of G protein signaling 7; COL22A1: Collagen Type XXII Alpha I Chain; MYC: V-Myc Avian Myelocytomatosis Viral Oncogene Homolog; MET: hepatocyte growth factor receptor; FGFR1: Fibroblast growth factor receptor 1; KIT: V-Kit Hardy-Zuckerman 4 feline sarcoma viral oncogene; EPHB4: EPH receptor B4; PIK3CA: Phosphatidylinositol-4,5-Bisphosphate 3-Kinase Catalytic Subunit Alpha; PTEN: Phosphatase and Tensin homolog; RON: Recepteur d’Origine Nantais or MST1R: Macrophage stimulating 1 receptor; ROS: reactive oxygen species; CTCs: circulating tumor cells; ECM: extracellular matrix; PI3K: phosphoinositide 3-kinase; ADAM-12: A disintegrin and metalloprotease-12; FAK: focal adhesion kinase; AKT: Protein kinase B; ERK: extracellular signal-regulated kinase; MSP: macrophage stimulating protein; MSP: macrophage stimulating protein; VCAM-1: vascular cell adhesion molecule 1 or CD106; VLA-4: Very Late Antigen 4 or integrin α4β1; CXCR4: C-X-C Motif Chemokine Receptor 4; CA: cellular automata; ODE: ordinary differential equations; PDE: partial differential equations; ABM: agent based models; BioFVM: Finite Volume Method for biological problems; RTK: receptor tyrosine kinase; HDAC: histone deacetylase; CAR-T: chimeric antigen receptor T.

CONFLICTS OF INTEREST

None.

FUNDING

This manuscript was supported by NIH/NCI grants P30CA033572 and 1U54CA209978-01A1.