Mathematical models of breast and ovarian cancers

Abstract

Women constitute the majority of the aging United States (US) population, and this has substantial implications on cancer population patterns and management practices. Breast cancer is the most common women's malignancy, while ovarian cancer is the most fatal gynecological malignancy in the US. In this review, we focus on these subsets of women's cancers, seen more commonly in postmenopausal and elderly women. In order to systematically investigate the complexity of cancer progression and response to treatment in breast and ovarian malignancies, we assert that integrated mathematical modeling frameworks viewed from a systems biology perspective are needed. Such integrated frameworks could offer innovative contributions to the clinical women's cancers community, as answers to clinical questions cannot always be reached with contemporary clinical and experimental tools. Here, we recapitulate clinically known data regarding the progression and treatment of the breast and ovarian cancers. We compare and contrast the two malignancies whenever possible in order to emphasize areas where substantial contributions could be made by clinically inspired and validated mathematical modeling. We show how current paradigms in the mathematical oncology community focusing on the two malignancies do not make comprehensive use of, nor substantially reflect existing clinical data, and we highlight the modeling areas in most critical need of clinical data integration. We emphasize that the primary goal of any mathematical study of women's cancers should be to address clinically relevant questions. WIREs Syst Biol Med 2016, 8:337-362. doi: 10.1002/wsbm.1343 For further resources related to this article, please visit the WIREs website.

Figure 1

“Breast Anatomy Female”: For the National Cancer Institute © 2011 Terese Winslow LLC, U.S. Govt. has certain rights.

Figure 2

“Female Reproductive System”: For the National Cancer Institute © 2009 Terese Winslow LLC, U.S. Govt. has certain rights.

Figure 3. The breast cancer growth models proposed by (A) Speer et al. and (B) Norton

(A) The first equation corresponds to classical Gompertzian growth kinetics, where N(t) is the number of tumor cells measure at time t after the start of tumor growth, A₀ is the initial growth rate, and α is the rate of growth decay. In this model, time is incremented in intervals of 5 days. A random number, r₁, between 0 and 1 is generated at each time interval, and compared to a predetermined value A₄, defined as the probability that α undergoes a change in a 5-day period. If r₁>A₄, the tumor continues growing at the previous rate. However if r₁<A₄, α is reduced by an amount depending on A₄, r₂, another randomly generated number between 0 and 1, and A₃, the predefined determinant of the amount of change in α. This process is illustrated in the second equation. Computations of the simulated growth curves are performed until either simulation time runs out, set at 40 years elapsed since the beginning of tumor growth, or until a lethal in silico tumor burden threshold is reached, set at N(t)=10¹² cells. The Speer et al. model beings with 1 cell at time 0 and uses predefined values of A₀, α, A₄, A₃. A₀ and α are expressed in units of days⁻¹, and A₃ and A₄ are dimensionless. Simulation time is measured in days. Each computed growth curve is generated using the same baseline parameter set. (B) The first equation corresponds to classical Gompertzian growth kinetics, where N(t) is the number of tumor cells measure at time t after the start of tumor growth set at t = 0, N(0) is the tumor starting size, b is the rate of growth decay, and N(∞) is the limiting size. To generate the probability distribution function of b, the proportion, P_L(t) of patients who have died by time t since the onset of symptoms after having reached the lethal tumor size, N_L, is generated from the 250 breast cancer survival curve data set in (Bloom et al.). Re-arranging the Gompertz, Equation (2) is obtained, where P_L(t_i) represents the proportion of the 250 cancers with growth decay rate b<b_i. The model is initialized with a set of initial values for N(0), N_L, and N(∞), and a least-squares based numerical algorithm is used to determine the mean and standard deviation of b. A randomly generated initial value for b_i is then chosen from the computed distribution to calculate t_i to ensure N(t_i) = N_L, and the process described in the second equation is repeated until the value of b_i that provides the best fit to P_L(t_i) is found.

Figure 3. The breast cancer growth models proposed by (A) Speer et al. and (B) Norton

(A) The first equation corresponds to classical Gompertzian growth kinetics, where N(t) is the number of tumor cells measure at time t after the start of tumor growth, A₀ is the initial growth rate, and α is the rate of growth decay. In this model, time is incremented in intervals of 5 days. A random number, r₁, between 0 and 1 is generated at each time interval, and compared to a predetermined value A₄, defined as the probability that α undergoes a change in a 5-day period. If r₁>A₄, the tumor continues growing at the previous rate. However if r₁<A₄, α is reduced by an amount depending on A₄, r₂, another randomly generated number between 0 and 1, and A₃, the predefined determinant of the amount of change in α. This process is illustrated in the second equation. Computations of the simulated growth curves are performed until either simulation time runs out, set at 40 years elapsed since the beginning of tumor growth, or until a lethal in silico tumor burden threshold is reached, set at N(t)=10¹² cells. The Speer et al. model beings with 1 cell at time 0 and uses predefined values of A₀, α, A₄, A₃. A₀ and α are expressed in units of days⁻¹, and A₃ and A₄ are dimensionless. Simulation time is measured in days. Each computed growth curve is generated using the same baseline parameter set. (B) The first equation corresponds to classical Gompertzian growth kinetics, where N(t) is the number of tumor cells measure at time t after the start of tumor growth set at t = 0, N(0) is the tumor starting size, b is the rate of growth decay, and N(∞) is the limiting size. To generate the probability distribution function of b, the proportion, P_L(t) of patients who have died by time t since the onset of symptoms after having reached the lethal tumor size, N_L, is generated from the 250 breast cancer survival curve data set in (Bloom et al.). Re-arranging the Gompertz, Equation (2) is obtained, where P_L(t_i) represents the proportion of the 250 cancers with growth decay rate b<b_i. The model is initialized with a set of initial values for N(0), N_L, and N(∞), and a least-squares based numerical algorithm is used to determine the mean and standard deviation of b. A randomly generated initial value for b_i is then chosen from the computed distribution to calculate t_i to ensure N(t_i) = N_L, and the process described in the second equation is repeated until the value of b_i that provides the best fit to P_L(t_i) is found.

Figure 4. The ODE-PDE framework used to model breast cancer development, treatment and recurrence, subsequently used to model radiotherapeutic strategies

The system of equations describing the interactions of tumor cells (n), extracellular matrix (f), and matrix-degrading enzymes (m) is illustrated in the first equation. Therein, the terms in the first equation correspond to cellular proliferation, random motility and haptotaxis, defined in the model as the movement of tumor cells according to gradients of chemicals in the tumor environment. The second equation corresponds to the extracellular matrix degradation by existing tumor cells. Lastly, the terms in the third equation correspond to the diffusion of matrix-degrading enzymes secreted by the tumor cells, the production of new enzymes and natural decay. The fourth equation represents the biologically effective dose , where n is the number of radiotherapeutic fractions administered, d is the dose delivered per fraction, α is the coefficient of single-hit DNA double-strand breaks, and β is the number of DNA single-strand break pairs that combine into forming double-strand breaks. d is measured in Gray. The surviving probability S, i.e. the proportion of cells that survive the radiation-induced damaged is modeled in the fifth equation.

Figure 5. (A) The four compartment cell cycle and resistance framework modeled by Roe-Dale et al.. (B) The two compartment cell cycle framework modeled by Panetta

(A) In the system of equations illustrated in (1), cells are separated by cell cycle status in two states, and drug sensitivity status in two other states, for a total of four possible states: either G1 or S sensitive (N₁) or resistant cells (N₃) and G2 or M sensitive (N₂) or resistant cells (N₄). In this model, resistant cells are defined as cells that express the activated MDR1 gene. The terms in each equation correspond to cell-specific constant transitions rates between the four compartments. The treatment equation for the Roe-Dale et al. model is illustrated by (2), where N represents the matrix-vector notation for the four different compartments whose temporal dynamics is modeled by (1), T_i is the corresponding treatment matrix for drug i, and m is the number of administered treatments with drug i at time intervals of τ hours. The fraction of cells surviving treatment with doxorubicin (T_A in the model) and with CMF (T_C in the model) are described in Equations (3) and (4), respectively. (B) In the system of equations illustrated in (1), P is the number of proliferating tumors cells, and Q is the number of quiescent tumor cells. Additional parameters include γ, the growth rate of proliferating cells, α, the transition rate from the proliferating to the quiescent compartment, δ, the natural proliferating cell death rate, β, the transition rate from the quiescent to the proliferating compartment, and λ, the natural quiescent cell death rate. All parameters in the model are assumed to be positive, and constant. Herein, the system outlined in (1) represents a linear system of ODEs modeling the dynamics of the proliferating and quiescent cell compartments. The function f(t) described in Equation (2) represents a step function describing the effects of the chemotherapeutic treatment, e.g. paclitaxel. The periodic function modeling the paclitaxel effects is assumed to target only the proliferating cell compartment. In its functional representation, s is the strength of the drug, a is the active drug time, T is the period of paclitaxel administration and n stands for the n^th administered drug dose.

Figure 5. (A) The four compartment cell cycle and resistance framework modeled by Roe-Dale et al.. (B) The two compartment cell cycle framework modeled by Panetta

(A) In the system of equations illustrated in (1), cells are separated by cell cycle status in two states, and drug sensitivity status in two other states, for a total of four possible states: either G1 or S sensitive (N₁) or resistant cells (N₃) and G2 or M sensitive (N₂) or resistant cells (N₄). In this model, resistant cells are defined as cells that express the activated MDR1 gene. The terms in each equation correspond to cell-specific constant transitions rates between the four compartments. The treatment equation for the Roe-Dale et al. model is illustrated by (2), where N represents the matrix-vector notation for the four different compartments whose temporal dynamics is modeled by (1), T_i is the corresponding treatment matrix for drug i, and m is the number of administered treatments with drug i at time intervals of τ hours. The fraction of cells surviving treatment with doxorubicin (T_A in the model) and with CMF (T_C in the model) are described in Equations (3) and (4), respectively. (B) In the system of equations illustrated in (1), P is the number of proliferating tumors cells, and Q is the number of quiescent tumor cells. Additional parameters include γ, the growth rate of proliferating cells, α, the transition rate from the proliferating to the quiescent compartment, δ, the natural proliferating cell death rate, β, the transition rate from the quiescent to the proliferating compartment, and λ, the natural quiescent cell death rate. All parameters in the model are assumed to be positive, and constant. Herein, the system outlined in (1) represents a linear system of ODEs modeling the dynamics of the proliferating and quiescent cell compartments. The function f(t) described in Equation (2) represents a step function describing the effects of the chemotherapeutic treatment, e.g. paclitaxel. The periodic function modeling the paclitaxel effects is assumed to target only the proliferating cell compartment. In its functional representation, s is the strength of the drug, a is the active drug time, T is the period of paclitaxel administration and n stands for the n^th administered drug dose.

Figure 6. The Liu el al. breast cancer stem cell model

The population dynamics of the cell types considered is illustrated in the above equations. Therein, x_i(t) is the number of cells measured at time t of type i, where i=0 is the cancer stem cell phenotype, i=1 is the progenitor cell phenotype, and i=2 is the terminally differentiated cell phenotype. p₀(p₁) is the probability that a cancer stem cell (progenitor cell) divides into 2 cancer stem cells (progenitor cells), q₀(q₁) is the probability that a cancer stem cell (progenitor cell) divides into 2 progenitor cells (terminally differentiated cells), v₀ and v₁ are the synthesis rates representing the transition rates from the cancer stem to the progenitor cell compartment and respectively, from the progenitor cell to the terminally differentiated cell compartment. Lastly, d_i, where i = 0, 1 or 2, represents the death rates of cells of type i.

Figure 7. The tumor growth model proposed by Brown and Palmer

Tumors in both early and advanced stages were assumed to grow exponentially. Specifically, the early stage tumor growth model is illustrated in Equation (1). Therein, a is the size at which a particular tumor is detectable by histopathology, b is the (exponential) growth rate constant and t₁ is the time since the tumor became detectable by histopathology. Detection thresholds for each individually simulated growth curve were set to match the corresponding value found in the collected tumor data set. The advanced stage tumor growth model is illustrated in Equation (2). Therein, c is the log value of the tumor size at disease progression from early to advanced stage (estimated from the Monte Carlo simulation of tumor life histories), d is the difference between the log values of the tumor size at empirical diagnosis obtained from the collected tumor data set and the log value of the size at progression from the generated simulation, and t₂ is the in silico measured time since progression.

Figure 8. The plasma biomarker temporal dynamics model of Hori and Gambhir

The change in the mass of the plasma biomarker with respect to time is equal to the difference between the influx of plasma biomarker shed by the tumor cells, u_T(t), healthy cells, u_H(t) and the outflux of biomarker from the plasma, q_EL(t), are as illustrated in the first equation. The rate of biomarker entry into the plasma is the sum of the input from tumor cells (as modeled in the second equation) and from healthy cells (as modeled in the third equation). Tumor cell growth is represented here by either the Gompertzian growth model (the fourth equation) or the exponential growth model (the fifth equation). The healthy cell population is assumed to remain constant throughout simulation time, and is set at N_H(t) = N_H,0.

Figure 9. The one-compartment ODE model of Kohandel et al.. Modeled are tumor growth, and surgical and chemotherapeutic treatments

Kohandel et al. considered one population of tumor cells, a non-cell cycle specific drug and various growth and cell-kill laws formulated in the following manner. The dynamics of the number of tumor cells at time t, N(t), is described by differential functional forms for the growth law, where f(N) is the tumor cell growth dynamics (e.g. f(N) = aN for the exponential growth law, where a is the constant proliferation rate), G(t,N) describes the effects of the drug on the system , and I(t = t_surgery) is an indicator function (equal to 1, if t = t_surgery, and 0 otherwise). Differential functional forms chosen for G(t,N) are provided in the second equation. Surgery is assumed to be instantaneous, and to remove a fixed fraction of exp(−k_s) of tumor cells, where k_s is the fraction of removed cells during surgery.