Optimization of additive chemotherapy combinations for an in vitro cell cycle model with constant drug exposures

Proliferation of an in vitro population of cancer cells is described by a linear cell cycle model with 𝑛 states, subject to provocation with 𝑚 chemotherapeutic compounds. Minimization of a linear combination of constant drug exposures is considered, with stability of the system used as a constraint to ensure a stable or shrinking cell population. The main result concerns the identification of redundant compounds, and an explicit solution formula for the case where all exposures are nonzero. The orthogonal case, where each drug acts on a single and different stage of the cell cycle, leads to a version of the classic inequality between the arithmetic and geometric means. Moreover, it is shown how the general case can be solved by converting it to the orthogonal case using a linear invertible transformation. The results are illustrated with two examples corresponding to combination treatment with two and three compounds, respectively.


Introduction
Recent decades have seen an increased interest in combination therapies to treat cancer [1]. Mathematical modeling of cancer growth and treatment plays an important role in the discovery and development of pharmaceutical compounds, as well as in increasing our understanding of biological processes [2]. Mathematical tools and prediction techniques are especially suitable to support the combinatorial explosion associated with testing different pairs, or even triplets, of anticancer drugs simultaneously and at different dose levels [3].
The mammalian cell cycle is central to cancer growth and therefore also to our understanding of how to treat cancer. The cell cycle is typically divided into four or five stages: 1 , , 2 , , and also 0 . Here, 1 and 2 are gap, or growth, phases that separate the stages , during which DNA synthesis takes places, and , which is when mitosis occurs. The stage 0 is a quiescent stage where a cell lays dormant in between cell cycles.
Chemotherapeutic compounds are often classified as cell cycle nonspecific, meaning that they target cells in all stages of the cell cycle, as well as quiescent cells, or cell cycle specific, meaning that they target cells in only one or a couple of the stages of the cell cycle [4]. Examples of cell cycle nonspecific drugs include platinum-based chemotherapies such as cisplatin and carboplatin. Examples of chemotherapeutic compounds that are cell cycle specific include enzymes such as asparaginase, which primarily target cells in the 1 phase, antimetabolites such as 5-fluorouracil and gemcitabine, which target cells in the phase, * Correspondence to: Fraunhofer-Chalmers Research Centre for Industrial Mathematics, Chalmers Science Park, SE-412 88 Gothenburg, Sweden.
topoisomerase including topotecan and irinotecan that target the 2 phase, and taxanes such as paciltaxel as well as vinca alkaloids such as vinorelbine that primarily target cells in the phase of the cell cycle [4].
Mathematical models that describe the growth of cancerous cells and tumors have been developed with varying complexity, taking into account many important biological features [5][6][7]. Perhaps the simplest model is described by a single difference or ordinary differential equation that captures growth that is either purely exponential, or which slows down or saturates as the tumor becomes large, e.g., in the Gompertz model [8,9]. More advanced models use systems of differential equations, and partial differential equations can be used to describe spatial or age-related aspects of tumor growth [10][11][12][13][14].
In this paper, we consider a simple mathematical model of the cell cycle model with stages. Early use of such a model to describe the growth of a population of cancer cells can be found in two papers by Takahashi [15,16]. Like many biological models based on ordinary differential equations, it can be recovered from a stochastic model where the time that a cell spends in each stage is assumed to be exponentially (or Erlang) distributed [17,18]. This assumption is consistent with some experimental results [19]. These models have been analyzed with respect to cell cycle kinetics without treatment [20][21][22][23][24], and have also been used to describe treatment with cell cycle specific drugs that target cells in certain stages of the cell cycle [25][26][27][28][29]. Similar models have also been used to describe combination treatments [30,31]. Cell cycle model for a population of cancer cells with states 1 and 2 , subject to combination therapy with two drugs, 1 and 2 , with concentrations 1 and 2 , respectively. Cells are transferred from state 1 with rate 1 and from 2 back to 1 with rate 2 during which mitosis occurs. Drug 1 acts cytotoxically on cells in state 1 with rate 11 and on cells in state 2 with rate 12 . Drug 2 acts cytotoxically on cells in state 1 with rate 12 and on cells in state 2 with rate 22 .
Our analysis of the linear cell cycle model incorporates combination treatment with an arbitrary number of chemotherapeutic drugs, whose drug actions may be different for cells in different stages of the cell cycle. In comparison, earlier works using a similar model have only considered special cases with one or multiple specific drugs [7,15,[30][31][32]. We first analyze the set exposure combinations that result in stability and refer to this as the shrinkage set. Then, we find an optimal treatment combination under the assumption of constant drug exposure by minimizing a weighted sum of the exposures subject to the constraint that the system is stable and therefore that the cancer cells will eventually be eradicated. The analysis rests on two main assumptions: (i) that cell killing is proportional to drug concentration, and (ii) that drug actions is additive, i.e., the model does not include general expressions for potential antagonistic or synergistic interactions.
Our results are illustrated using two examples. The first example combines the drugs 5-fluorouracil and vinorelbine, which act primarily on the and phases of the cell cycle, respectively. The second example incorporates a third compounds, irinotecan, which acts mainly on the 2 phase of the cell cycle, and considers the problem of finding optimal triple combinations that result in a stable or regressing population of tumor cells.
We end the introduction with an analysis of the special case with a cell cycle model with two stages, subject to combination treatment with two chemotherapeutic compounds. Consider a simple model of the cell cycle with two states 1 and 2 . Cells in 1 travel to 2 with rate 1 . At 2 mitosis occurs and the two daughter cells are transferred back to 1 with rate 2 where the cycle repeats. Two chemotherapeutic drugs, 1 and 2 , induce cell death via apoptosis depending on the concentration of the drugs, denoted by 1 and 2 , respectively. Assume that drug action is linear and that apoptosis is induced for cells in state , due to drug , with rate ≥ 0. The model is illustrated in Fig. 1. Growth of the cell population subject to chemotherapeutic provocation with the two drugs is then described by the following system of differential equationṡ 1 = −( 1 + 11 1 + 12 2 ) 1 + 2 2 2 , 1 (0) = 01 , 2 = 1 1 − ( 2 + 21 1 + 22 2 ) 2 , where 01 and 02 are the number of cells in states 1 and 2 at time = 0, respectively. Next, we investigate the stability of the system (1) to determine which combinations ( 1 , 2 ) lead to eradication of the cancer cell population. The system (1) has the system matrix ] .
(2) It is well-known that, for a 2 × 2-matrix, the system is stable if and only if tr < 0 and det > 0. We note that the condition on the trace always holds, and that the determinant is given by The hyperbola intersects the coordinate axes in two points, ( * 1 , 0) and (0, * 2 ), corresponding to the minimum concentration of either compound that, when given as monotherapy, results in stability. The determinant condition, det = 0, for monotherapy becomes det = 1 2 2 + ( 1 2 + 2 1 ) − 1 2 , = 1, 2, which is a second-order equation with exactly one positive root, given by * = −( 1 2 + 2 1 ) + √ ( 1 2 + 2 1 ) 2 + 4 1 2 1 2 so that any > * leads to eradication. In general, we define the set which we refer to as the shrinkage set. The boundary, , separates the first quadrant of the 1 2 -plane into a region corresponding to stability and a shrinking cell population, and a region corresponding to instability and population growth. As noted above, is (part of) a hyperbola, and the set is therefore convex. An illustration of this set is shown in Fig. 2.
We consider a linear objective function ( 1 , 2 ) ∶= 1 1 + 2 2 , where the coefficients correspond to some cost, e.g., related to toxicity and side-effects, associated with the drugs, respectively. This leads to the convex optimization problem which can be solved using the standard technique with Lagrange multipliers.

Methods
Note: Throughout the rest of this paper, inequalities for vectors and matrices are to be interpreted element-wise, e.g., ≥ for matrices of equal dimensions, or ≤ for vectors of equal length. We describe in vitro growth of a population of cancer cells with a general linear cyclic model with states as a description of the cell cyclė where represents the number of cells in the :th stage, 0 is the initial number of cells in stage , and are the transfer rates between stages. At the last stage, , mitosis occurs and the two daughter cells are placed in state 1 and re-enter the cell cycle. Unlike normal cells, we assume that the cancerous cells grow unimpeded, and never enter a quiescent stage 0 to wait for external signals to commence proliferation. The system (9) is linear and autonomous and can be written in matrix forṁ where = ( 1 , … , ) is the vector of states, 0 = ( 01 , … , 0 ) is the initial vector, and the matrix is given by It is easy to see that if 0 ≥ 0 the solution ( ) remains in the nonnegative orthant, , for all ≥ 0. Take ∈ . Then from (9) we have thaṫ≥ 0 for all such that = 0. Hence, the vector field never points outside and ( ) is trapped inside.
We extend the cell cycle model (10) to account for multiple drug provocations, by introducing a control vector = ( 1 , … , ) where represents the concentration of drug , assumed to be non-negative and constant over time, which is plausible for an in vitro setting. We assume that drug action is linear and may be different for different stages in the cell cycle. Although drug action is in general sigmoidal, a linear approximation can often be justified for a range of exposure levels such that saturation effects are sufficiently small [33]. Moreover, we assume additive effects for combinations of drugs that target the same stage of the cell cycle. This is a natural assumption given that we consider a general case and not specific compounds with known mechanisms. For each , let be a non-negative diagonal × -matrix, whose entries , represent the efficiency of drug at killing cells in state . This model is illustrated in Fig. 3. The model is described by the following system of differential equationṡ Eq. (12) is a linear system with respect to both states and drug concen- Definition 1 (Drug Action Matrix). We define the drug action matrix to be the matrix whose entries are given by = , .
The entries correspond to the drug action of drug on state . Using the drug action matrix, can be written We investigate the asymptotic behavior of the solution ( ) to (12) by considering the eigenvalues of . In particular, we are interested in the case when all eigenvalues have negative real part, which implies that tumor will shrink over time and eventually be eradicated. First, we note that can be written = − , for some ≥ 0, where is a non-negative irreducible matrix. It follows from the classical Perron-Frobenius theorem (see, e.g., [34]) that has a simple largest real eigenvalue , with a corresponding eigenvector with positive components. Thus, the asymptotic growth rate will be determined by and the corresponding distribution of cells into the different states will be determined by . Our objective is therefore to find those such that becomes negative. We make the following definition Definition 2 (Shrinkage Set). For any square matrix whose entries depend on a control vector , define the shrinkage set as the set of all non-negative such that all eigenvalues of have negative real part.
Before analyzing the shrinkage set for we make the following observation.
Gerschgorin's circle theorem states that the eigenvalues of a matrix are contained inside the Gerschgorin discs in the complex plane [35]. Applying the theorem to the columns of gives the discs (− − ∑ , ), = 1, … , − 1 and (− − ∑ , 2 ), where the first and second argument denotes the origin and radius of the disc, respectively. Only the last of these discs can intersect the positive half of the complex plane. Moreover, increasing the values of does not change the radii of the Gerschgorin discs, but only pushes them further into the left half of the complex plane. If ∑ ≥ all discs are located in the left half of the complex plane, i.e., no eigenvalue of has a positive real part. Applying the theorem to the rows of gives the similar condition ∑ ≥ 2 − 1 . Note also that if there exists such that = 0 for all , i.e., there exists a state that is not acted upon by any drug, then even if all other discs are pushed far into the left half of the complex plane, the corresponding Gerschgorin disc does not change and therefore provides a bound on how quickly the tumor can shrink.
The matrix has non-negative off-diagonal entries and is thus a Metzler matrix (c.f. [36]). The additive inverse of such a matrix is known as a Z-matrix. Z-matrices have been studied extensively, as they can be written − , where is a scalar, and is a non-negative matrix. Metzler matrices and Z-matrices occur frequently in applications, e.g., in the study of linear compartment models [37]. What follows is a presentation of selected parts of the theory of Z-matrices, which are useful in establishing our main result.
Recall that a matrix is said to be positive stable if all of its eigenvalues have positive real part. A positive stable Z-matrix is known as an M-matrix. Hence, the question of stability for a Metzler matrix such as in Eq. (13) can be recast as a question of whether its additive inverse is an M-matrix. In 1979, Berman and Plemmons published a collection of a large number of properties that are equivalent to a Z-matrix also being an M-matrix [38]. Some of these properties are useful to us and are summarized in the following theorem Theorem 1 (Berman and Plemmons, 1979). Let ∈ R × be a Z-matrix. The following are equivalent to also being an M-matrix Recall that a principal minor of a square matrix is the determinant of a submatrix obtained by removing a number of rows and the corresponding columns. The leading principal minors are the determinants of the submatrices obtained by successively removing rows and column, starting with the last row and column. Note that the determinant of a matrix counts as a principal minor to itself. The sum of two Z-matrices is clearly a Z-matrix. However, it is well-known that the sum of two M-matrices need not be an M-matrix. Some results regarding convex combinations of M-matrices have been presented, c.f., Fan [39], Horn and Johnson [34], and Stipanović and Šiljak [40]. We shall use a result that follows immediately from a theorem originally proved by Cohen [41] although the version presented here is due to Friedland [42] Theorem 2 (Cohen, 1978). Let 1 and 2 be two M-matrices such that 1 − 2 is a diagonal matrix. Then, for any ∈ [0, 1], the matrix Note that this result follows from Theorem 1 when 1 − 2 ≥ 0. An immediate consequence is the following proposition Note that the proposition gives convexity for the shrinkage set of any compartment model where the control vector stimulates or induces output from the system. Note also that convexity is preserved if we replace the linear functions ( ) with any concave functions of , e.g., saturable functions of Michaelis-Menten type [43].
Suppose now that we want the tumor described by (12) to be eradicated, i.e., the matrix should be stable, while at the same time minimizing the metabolic strain induced by the drugs. A simple and common way to express this is based on the (weighted) total drug exposure ( ) given by where ≥ 0 is a vector of weights that reflect the relative toxicity of the drugs . Such objective functions have been used in various optimal control problems of cancer therapies, c.f. the book by Schättler and Ledzewicz [44].

Results
We begin by presenting the main results, with proofs, followed by two illustrative examples.

Lemma 1. The shrinkage set for the matrix given by Eq. (13) is the convex set given by
Proof. We first note that convexity follows from Proposition 1. To find a formula for the shrinkage set , it follows from Theorem 1 that we only need to study when the leading principal minors of − are positive. All leading principal minors except det(− ) are trivially positive since they are determinants of positive diagonal matrices. Therefore, is determined by the condition det(− ) > 0. An application of Leibniz's formula for determinants gives condition (16). □ Note that can also be expressed using the drug action matrix The Eq. (17) can be given a probabilistic interpretation. First, recall that a system such as (12) can be derived from a probabilistic model where the transfer time, , from state to the next state is an exponentially distributed random variable with parameter . Similarly, cell death in state due to drug can be assumed to be exponentially distributed with parameter . By independence, the time until cell death for a cell in state due to any of the drugs is exponentially distributed with parameter ∑ . We denote this time by . A cell in state survives until the next state if < . The probability of this event is given by Letting be a random variable such that = 1 if a cell survives the entire cell cycle, and = 0 if the cell dies somewhere along the way. Then, by independence, the probability of surviving is given by Since each completed cell cycle results in two new cells, the tumor is expected to shrink if P( ) < 1∕2, which gives the same condition as Eq. (17). We first use this lemma to prove the following intuitively obvious fact about monotherapy Proposition 2 (Monotherapy). Consider the monotherapy case with only one compound, whose tumor concentration we denote by , acing on each of the states 1 , … , according to the diagonal matrix . Then, there exists exactly one value * such that any > * ensures stability of .
Proof. Define the polynomial given by and note that (0) < 0. By Lemma 1, is stable precisely for those for which ( ) > 0. Rewriting on the form we note that 0 < 0 and > 0 for all ≥ 1. It follows from Descartes's rule of sign that has exactly one positive root, * , and consequently ( ) > 0 for all > * . □ Proposition 1 expresses the fact that for a single compound, there exists exactly one concentration * such that any concentration below it will lead to tumor growth, whereas any concentration above will lead to tumor shrinkage. Finding an explicit formula for * is non-trivial, but it can be expressed using the trick from the proof of part (iii) of the next proposition involving combination therapy.
The next proposition considers the combination therapy case, with drugs. We characterize the optimal solution, given a linear objective function in and using the stability condition as a constraint. We optimize over the closed set̄, which will give an optimum * ∈ . Then, for any vector > 0 we have that * + ∈ , which will ensure stability. We also assume, without loss of generality, that the drug action matrix does not contain rows of all zeros, corresponding to states without drug action, since it is clear from Eq. (17) that the only influence from such states comes from the corresponding , which can be simplified from Eq. (17).

Proposition 3 (Combination Therapy). Let the shrinkage set be defined by Eq. (17) and consider the optimization problem
where = ∕ , and denotes the geometric mean. The optimum is then given by * = and the optimal value becomes

(iii) If has full rank, the problem can be converted to case (ii) by a linear invertible transformation.
Proof. We first prove part (i). Note that, due to Lemma 1, the domain can be written Since the objective function is linear, any optimum * must be located on the boundary ( ) = 0. If is rank-deficient with rank − for some positive integer , then for each ∈ ( ) the solution set to = is a -dimensional affine space. If * > 0 is an optimum, then the -dimensional tangent spaces of and at * must coincide and, since is linear, be a level set to . Now, since = * is a -dimensional affine space, there must exist a point̃, with zeros, in this level set, that also satisfies (̃) = 0. This proves (i).
To prove (ii), note that when is diagonal, the change of variables = + gives rise to the equivalent optimization problem, which also occurs, e.g., in one of the standard proofs of the inequality between the arithmetic and geometric means, and ≥ , where = ∕ . This is a convex optimization problem, hence any local minimum must also be a global minimum (see, e.g., [45]). Since the objective function is linear and the optimum therefore is attained on the boundary, we can assume that Eq. (28) holds with equality. Define the Lagrangian The Karush-Kuhn-Tucker conditions, which are necessary for optimality by the linear independent constraints qualification (see, e.g., [45]), are given by ≥ 0.
Let = ∏ =1 . Multiplying Eq. (31) with and using Eq. (32) gives Moreover, taking the product of Eq. (31) over all and using Eq. (32) gives Combining Eqs. (35) and (36) gives the solution * = Finally, if * > 0, then = 0. Letting = ∕ gives * = The condition (23) follows by noting that it precisely the condition that ensures positivity in Eq. (38). This proves (ii). Finally, we prove (iii) by constructing the desired linear invertible transformation. Note that the claim is obvious when is invertible, since the change of variables = gives case (ii) with new cost vector = − . If has full rank, but does not have enough columns to span R , can be augmented with Euclidean unit vectors +1 , … , such that̃∶ = [ , +1 , … , ] is invertible. Letting the additional columns of̃correspond to control variables +1 , … , with costs +1 , … reduces the problem to the invertible case with the new variables = and costs =̃− , which gives the analog to Eq. (37) * = from which it follows Finally, we note that by choosing the fictitious costs such that they satisfy the equations where =̃− is the transformed cost vector, ensures that * = 0 for = + 1, … , so that the fictitious control variables do not appear in the solution. □ Case (i) states that if we have more compounds than targets, then we do not lose anything by getting rid of some of the compounds. Case (ii) presents the diagonal case where each compound has a unique single target. This could also be called the orthogonal case, since it covers any orthogonal matrix up to scaling and a reordering of the compounds. An important insight expressed in case (ii) is that there is intrinsic benefit from combining compounds with different targets, which is precisely due to the convexity of the shrinkage set established in Lemma 1 and Proposition 1. Case (iii) relates the general case back to case (ii) via a linear bijection.
We have the following obvious corollary to Proposition 3, case (ii), which can also be applied to case (iii), after a linear change of variables. Corollary 1. If * = 0 for some , then equation (24) simplifies to * = Proof. The proof is a direct calculation, letting * = 0 in Eq. (37) and solving for , then plugging this value back into Eq. (38) for ≠ . □ We can also consider, as a special case, the situation where all compounds are 'equal' in the sense that the are equal. This gives an estimate of the inherent benefits of targeting different states.
Proof. The proof is a direct calculation. The limit follows, e.g., by an application of l'Hopital's rule. □ For a single drug, = 1, Corollary 2 gives the objective value ( * ) = , whereas in the limit as → +∞ we have ( * ) = ln 2 ≈ 0.7 . Thus, Corollary 2 states that combination therapy with additive drugs that target different stages of the cell cycle can at most reduce the cost compared to single-agent treatment by 30%.

An example with a two-drug combination
We consider a cell cycle model with four states, 1 , 2 , 3 , and 4 , corresponding to the stages 1 , , 2 , and of the cell cycle, and transfer rates 1 , 2 , 3 , and 4 . In order to obtain biologically plausible parameter values, we assume phase durations of 12, 6, 6, and 1∕2 hours for the phases 1 , , 2 , and , respectively [46]. Since the average time in the i:th phase is given by 1∕ , we therefore set 1 = 1∕12, 2 = 1∕6, 3 = 1∕6, and 4 = 2 h −1 . With these parameter values, the largest and only positive eigenvalue of the corresponding system matrix is given by ≈ 0.032. The corresponding normalized eigenvector is ≈ (0.55, 0.23, 0.20, 0.02). Thus, the initial proportion of cells in each phase is assumed to be given by . We now consider treatment with two cell cycle specific compounds that target different phases of the cell cycle. 5-fluorouracil is an antimetabolic agent that interferes with DNA and RNA synthesis, and therefore acts mainly on the phase of the cell cycle, that is used to treat a variety of cancers including stomach cancer, colon cancer, pancreatic cancer and breast cancer [4]. Vinorelbine is a vinca alkaloid that interferes with microtubule assembly and therefore acts primarily on the stage of the cell cycle [4]. Vinorelbine is used to treat different cancers including breast cancer and non-small cell lung cancer [47].
Let denote the drug action parameter of 5-fluorouracil on cells in phase, and let denote the drug action parameter of vinorelbine on cell in phase. Moreover, let be the cost, or toxicity, or 5fluorouracil and be the toxicity of vinorelbine. Combinations of these compounds satisfy case (ii) of Proposition 3, i.e., the orthogonal case. Monotherapy with 5-fluorouracil requires exposure levels above 2 ∕ to achieve stability and a shrinking population of cancer cells, whereas the equivalent quantity for vinorelbine monotherapy is 4 ∕ . Note that since the phase is much shorter than the phase of the cell cycle, the potency parameter would need to be much larger than in order to achieve stability at similar exposure levels.
To compute biologically reasonable values for and we use reported in vitro geometric mean IC50 values across a large number of cell lines of 34.8 μM and 0.0217 μM for 5-fluorouracile and vinorelbine [48], and then use our model to find the corresponding and such that the cell population after 72 h of single-agent treatment is exactly half of the untreated cell population at the same time point. This yields the values = 19.4 μM −1 h −1 and = 0.00115 μM −1 h −1 . The shrinkage set, defined in Lemma 1, is given by This set is shown in Fig. 4. In terms of Eq. (23) we have that 5 = 2 ∕ and = 4 ∕ , and therefore ( ) = √ 2 4 ∕ √ . It will be optimal to use only 5-fluorouracil if 5 ≥ 2 and it will be optimal to use only vinorelbine if ≥ 2 5 . Therefore if or equivalently it will be optimal to use both compounds and the optimum will be given by * 5 and * = 4 After inserting the values for 2 , 4 , , and , Eq. (46) defines a condition on the costs and for which it will be optimal to use both compounds. This set in shown in Fig. 5. The light blue region shows the cost pairs ( , ) such that it will be optimal to use both compounds, whereas the white regions on either side indicate for which cost pairs it will be optimal to use only 5-fluorouracil or vinorelbine.

An example with a three-drug combination
We consider the same two compounds as in the previous example, 5-fluorouracil and vinorelbine, but this time also include irinotecan. Irinotecan is a topoisomerase I inhibitor, which stops the cell at the G2-checkpoint which eventually triggers apoptosis. Irinotecan is used to treat primarily colorectal cancers [49]. We assume that irinotecan acts on cell in the G2 phase with rate , and let its cost be denoted by . The cell cycle model subject to triple combination treatment with 5fluorouracil, vinorelbine, and irinotecan is shown in Fig. 6. To compute we take the same approach as before and use the model to find , , respectively. It is assumed that 5-fluorouracil acts on cells in the phase with parameter , vinorelbine acts on cells in phase with parameter , and irinotecan acts on cells in state 2 with rate parameter .
such that the cell population after 72 h of irinotecan treatment at reported IC50 concentrations is exactly half the size of the untreated cell population. The reported geometric mean IC50 for irinotecan across a large number of cell lines is 13.8 μM [48], which corresponds to = 0.0029 μM −1 h −1 .
The shrinkage set, defined in Lemma 1, for combinations of 5fluorouracil, vinorelbine, and irinotecan is given by the condition and is shown in Fig. 7. Now let = 3 ∕ , and ( ) = ( 5 ) 1 3 . Then the conditions that ensure that it is optimal to use all three compounds are given by Using the estimated values for 1 , 2 , 3 as well as , , this gives rise to a region in cost space consisting of triples ( , , ) for which it is optimal to use all three compounds. This region is shown in Fig. 8. It is clear from the figure that in order for the optimum to use all three compounds we would need ∼ 10 2 and in turn ∼ 10 2 .

Discussion
We analyze stability conditions for a linear cell cycle model subject to chemotherapeutic combination treatment. Such analyses are common for single-agent treatments and lead to the definition of threshold concentrations such that any concentration above the threshold results in a shrinking tumor or population of cancer cells [50] and have also been derived for combination treatments [32,51]. We prove the existence of such a concentration for our model in Proposition 2. However, finding an explicit formula for the general case where the single-agent treatment may act differently on different phases of the cell cycle is difficult. One way to approach this problem is to apply Proposition 3(iii), and the proof thereof, when there are more targets than compounds. For combination treatments, stability analysis leads to the definition of a shrinkage set consisting of all exposure combination that result in a shrinking cancer cell population. The geometry of the shrinkage set is interesting since it is related to the inherent benefits of combination therapy [51]. For the cell cycle model in this paper, where the anticancer agents are assumed to be additive, we show in Lemma 1 and Proposition 1 that the shrinkage set is always convex, and that this will be true for general compartment models subject to treatment with cell killing drugs.
We use the shrinkage set as a constraint to optimize combination treatments under the assumption of constant drug exposure. Constant exposure is achievable in an in vitro setting which is our primary consideration in this paper, but it is less feasible in vivo. Instead, constant exposure should be viewed as an approximation that is suitable for drugs that can be administered via continuous infusion, or which have frequent dosing and/or slow enough clearance such that exposure levels can be maintained at fairly stable levels. Our main result, given in Proposition 3, provides a formula for the optimal exposure combination as well as conditions for when one or several compounds are superfluous. We believe that these results can be useful when considering new combinations treatments with compounds for which cell cycle specificity is at least somewhat understood. Moreover, we believe that Corollary 2, which provides a bound of at most 30% cost reduction compared to single-agent treatment, to be a relevant result that can help determine whether additive combination treatments are feasible or whether synergistic compounds with drug interactions are necessary to achieve a desired cost reduction.
Optimization and optimal control in the more general case of cancer and cell cycle models with time-varying exposure have been analyzed in numerous papers including the works of Kimmel, Swierniak, Ledzewicz, and Schättler to name a few [44,[52][53][54][55][56]. The objective to be minimized is typically a weighted 1 -or 2 -norm of the drug concentrations and tumor volume, combined with an endpoint penalty at a final time, which can be either fixed or free. Compared to the case with constant exposure, which is appropriate for the in vitro setting, optimal control problems are more suitable for the clinical setting [57]. Moreover, although optimal controls can be computed numerically for specific choices of parameter values, it is often very difficult to obtain explicit analytical solutions. One important result, which has been proved to hold for a general case of linear cancer models, states that optimal controls are always ''bang-bang'', i.e., at any time point the exposure level should be set either to a maximum allowed value, or to zero [44]. However, a complete determination of the switching times, when the exposure should be changed from its maximum value to zero or vice versa, is highly nontrivial.
The optimization problem with constant drug exposures that we consider in this paper is much simpler than the related optimal control problems and can therefore be analyzed analytically. Another difference is that whereas optimal control problems consider the full dynamics of the model as a constraint, we only consider the sign of the eigenvalues and therefore mainly consider the asymptotic behavior of the system. One important result by Clairambault et al. states that a linear model given by a Metzler matrix, with time-varying periodic coefficients, has worse stability properties than the corresponding system with time-averaged constant coefficients [58]. This provides justification for our restriction to only consider constant drug exposures.
The model considered in this paper describes an in vitro population of cancer cells. In order to describe in vivo tumor growth, the model would need to be modified to account for spatial heterogeneity in particular. This can be done using partial differential equations, or an ordinary differential equation approach such as the one by Evans et al. [59] who describe a tumor consisting of an outer shell and an inner core with the possibility to transfer cell mass between the two. See the paper by Checkley et al. for an example where the shell model was integrated into a cell cycle model [60].

Conclusions
This paper analyzes a simple in vitro model of a population of cancer cells that captures two important aspects of cancer treatment: (i) the cell cycle, and (ii) combination treatment with two or more chemotherapeutic compounds. We find an explicit condition for stability of the model, and show that the corresponding shrinkage set is convex, which is useful for understanding which exposure combinations may produce a stable or shrinking cancer cell population.
Moreover, our main result, Proposition 3, provides a condition for when compounds are redundant, and a solution formula for the case where all compounds are used. These results could potentially be useful for practical applications such as selecting compounds for additional experiments, as well as for more mathematical endeavors such as sensitivity analyses and to understand for which sets of parameter values a particular combination is optimal.
Finally, since our analysis is performed for general additive combination treatments, it has the potential to be used and reused to analyze many different chemotherapeutic scenarios. It could also be used as a baseline case for additivity, which may then help determine if a particular combination is synergistic or antagonistic.