 Research
 Open Access
 Published:
A patientspecific therapeutic approach for tumour cell population extinction and drug toxicity reduction using control systemsbased doseprofile design
Theoretical Biology and Medical Modellingvolume 10, Article number: 68 (2013)
Abstract
Background
When antitumour therapy is administered to a tumourhost environment, an asymptotic tapering extremity of the tumour cell distribution is noticed. This extremity harbors a small number of residual tumour cells that later lead to secondary malignances. Thus, a method is needed that would enable the malignant population to be completely eliminated within a desired timeframe, negating the possibility of recurrence and druginduced toxicity.
Methods
In this study, we delineate a computational procedure using the inverse inputreconstruction approach to calculate the unknown drug stimulus input, when one desires a known output tissueresponse (full tumour cell elimination, no excess toxicity). The asymptotic extremity is taken care of using a bias shift of tumourcell distribution and guided control of drug administration, with toxicity limits enforced, during mutuallysynchronized chemotherapy (as Temozolomide) and immunotherapy (Interleukin2 and Cytotoxic Tlymphocyte).
Results
Quantitative modeling is done using representative characteristics of rapidly and slowlygrowing tumours. Both were fully eliminated within 2 months with checks for recurrence and toxicity over a twoyear timeline. The dosetime profile of the therapeutic agents has similar features across tumours: biphasic (lymphocytes), monophasic (chemotherapy) and stationary (interleukin), with terminal pulses of the three agents together ensuring elimination of all malignant cells. The model is then justified with clinical case studies and animal models of different neurooncological tumours like glioma, meningioma and glioblastoma.
Conclusion
The conflicting oncological objectives of tumourcell extinction and host protection can be simultaneously accommodated using the techniques of drug input reconstruction by enforcing a bias shift and guided control over the drug dosetime profile. For translational applicability, the procedure can be adapted to accommodate varying patient parameters, and for corrective clinical monitoring, to implement full tumour extinction, while maintaining the health profile of the patient.
Background
Though there has been considerable efforts in exploring newer modalities of cancer treatment in the last several decades, hopes have been belied by a fundamental reason, though not fully appreciated, namely the difficulty of eradicating a tumour mass due to the nature of reaction kinetics that govern the interaction of tumour cells with the therapeutic agents administered. These interactions are governed by firstorder reaction processes (chemotherapy) and enzyme saturation kinetics (immunotherapy) [1]. Both these prescribed treatments cause an exponential decay of the tumour cell population leaving a finite definitive number of tumour cells at the asymptotic extremity of tumour versus drug distribution curve. By its very nature, an asymptotic tail implies the tumour cell population curve will contact the horizontal axis and become zero only when either the time duration or the drug dose is infinite. For instance, a typical tumour may have 10^{10} cells, and the surviving fraction, SF, of tumour cells after administration of a drug concentration D is SF = exp (−α D). Using the standard values of α = 0.02 and D = 150 mg for bleomycin (maximum dose tolerated) [2], the equation yields 1,765 cells surviving. This elucidates why many tumours recur, after appearing to initially shrink or regress under therapy.
An imminent question in neurooncology is the efficacy of the drug to be able to cross the “blood–brain barrier” (BBB) to enter the brain. A recent progress is the development of DNA alkylating drugs such as Temozolomide (TMZ) which while being highly targeted, can effectively cross the BBB [3]. However TMZ is known to be contraindicated in some cases such as in patients with severe myelosuppression. A few other common side effects are nausea and vomiting which are selflimiting or readily controlled with standard therapy. Temozolomide is a teratogenic compound and thus should not be used during pregnancy. It might very rarely leads to acute respiratory failure. However, no drugrelated adverse CNS effects or alopecia are known to occur with temozolomide [3].
To complement the effects of chemotherapy, immunotherapy, is also being considered for its synergistic carcinolytic effects. Recent incisive findings of Wheeler et al. [4] indicate that a combinatorial therapy design utilizing immunotherapy with chemotherapy reduces tumour volume by 50% and appreciably extends the average 1 year survival duration of glioblastoma patients, which cannot be done by either chemotherapy or immunotherapy alone. An important aspect of immunotherapy is to utilize cytotoxic Tlymphocytes (CTL) which are CD8 + T cells, that are known to be carcinolytic. Activated CTL can be generated by administering immunomodulating factors,. This can be done by two means: (i) administering cytokines such as Interleukin2 (IL2) which can cross the blood–brain barrier, and (ii) injecting cellular agents such as tumourinfiltrating lymphocytes (TIL) prepared beforehand by sensitizing Tcells of the patient’s blood against the tumour (biopsy tissue). These TILs proceed satisfactorily to the tumour mass in the brain parenchyma. Also, IL2 is well known to stimulate recruitment and proliferation of cytotoxic Tlymphocytes, such as CD8+ T cells, suggesting a novel neurooncological approach [5].
A critical aspect of anticancer therapy is maintaining protection of the healthy (nonmalignant) cells of the affected as well as neighboring unaffected tissue within tolerable limits. In cases of toxicity due to high doses of chemotherapeutic or immunotherapeutic agents, the patient may succumb to opportunistic pathogenic infection occurring due to the massive reduction of the normal protective immune cells by the agents, thus leading to septicemia, multiorgan failure and eventually death. The normal cell protection can be characterized by the concentration of circulating lymphocyte population in the blood [6]. These lymphocytes generate antibodies against infection, activate Bcells, Tcells and Natural Killer (NK) cells against pathogens. For ensuring that normal tissue protection is maintained during therapy, we have applied constraints or bounds in the feedback function of the drug dosages in the therapy control system, so that there is always a selfregulated judicious combination of the doseandtime schedule of the antitumour agents, such that the circulating lymphocytes do not fall below a tolerable threshold value. The therapeutic agents have compensatory effects on the lymphocytes, as the chemotherapy input damages the lymphocytes, while the immunotherapy drug protects the lymphocytes. A measure of the tissue damage inflicted by antitumour agents can be described by the toxicity cost function J, wherein cytotoxicity, a secondorder effect [7], depends on secondpower of drugdosage:
Here u_{ 1 }, u_{ 2 }, …. are the levels of the different antitumour agents, while B_{ 1 }, B_{ 2 }, …. are the weighting factors of the different agents. We use this principle to suitably orchestrate the temporal schedule of the drugs, so that that toxicity is minimized.
We may mention that various attempts at modeling the immune system interaction with neoplastic tumours have been previously made [8–10]. These models have efficiently characterized the computational dynamics of drug versus tumour interaction via the immune system. Using the background of the existing models, in our model we have tried to delineate the kinetics and dynamics of immune modulation responsible for the paradoxical clinical phenomenon of tumour dormancy, prolonged arrest and oscillations of tumoursize [11]. A unitary approach to the dual behaviour of tumour progression and tumour regression has recently been explained [12], where the neoplastic process has been elucidated as systems biologybased abnormality. The tumour regression approach that we report in the present work is to our knowledge, the first endeavor to elucidate a quantitative methodology to delineate the dosetime profile of administration of the antitumour agents (chemotherapy, interleukin, lymphocytes) with neuroncological cases as examples, so as to enforce the tumour cell population to zero, thus enabling full tumour elimination. For this, we develop an interdisciplinary approach, utilizing input reconstruction analysis and bias shift.
Methods
Inverse construction of drug input for obtaining desired tumour response
In conventional quantitative mathematical models, the inputs in terms of therapeutic agent concentration are substituted in the model (differential equations of tumour cell population) to solve for tumour cell population at different dosages and time durations. This is a “forward direct solution”, whereby, given the stimulus or input (concentrations of the N number of drugs, D_{ N }), one calculates the response or output (remaining tumour cell population T), by considering the successive intermediary stages of celldrug interactions (X, Y, Z):
Typically, in eq. (1), the investigator starts with varying dosages of the drug combinations D_{ N }, and obtains the resulting values of the surviving tumour cell population T as one changes the dosage profiles of the drugs. Using this information, the specific dosing schedule of the drugs, D_{ N }^{†} is selected, which would lead to minimization of the number of remaining cancer cells, T_{ min }. The sites X, Y and Z are referred respectively as intrinsic (input), intermediate and extrinsic (output) compartments (Figure 1). The extrinsic compartment Z is related to the external output of the therapy system (i.e. remnant tumour cell population T). Likewise, the intrinsic compartment X relates to the internal control input coming from the therapy system, i.e., administration of therapeutic agents, D_{ N }.
To circumvent the problem of T_{ min } i.e. the residual cancer cells that may cause tumour recurrence, we suggest the “backward inverse solution”, whereby one starts first with the desired objective, i.e. response or output (viz., tumour cell population becoming zero in a desired finite time period, t_{ 1 }), and, then work backwards through the sequential stages inversely, namely Z, Y, X to derive the drug concentrations to attain the above mentioned objectives:
So, in the equation (2), one puts T = 0 at time t_{ 1 }, which is used to calculate the temporal profile of dosages of the combination drugs (D_{N}^{†}). Thus, if one administers the temporal dosages profile (D_{N}^{†}), then the tumour cell population definitively becomes zero in time t_{ 1 }, indicating tumour elimination. The perspective of inverse analysis approach [eq. (2)] has been well used in other fields in the form of stimulus reconstruction approach. This enables one to accurately quantify an antecedent unknown stimulus, using the knowledge of the observed response pattern. In our case, the methodology enables one to reconstruct the input drug stimulus, which will accurately produce complete tumour elimination as the output in finite time.
Problem of asymptotic drug reaction kinetics
The celldrug interactions are governed by firstorder chemical reaction rate kinetics [1] using which the tumour cell population T at time t decreases exponentially in a monotonic manner:
where T_{ 0 } is initial population, T′ is the time derivative of T, and κ_{ T } is the drug’s cytotoxic rate constant (Figure 2A). As explained before, there persists asymptotically a residual tumour cell population, which now needs to be dealt with. Hence, one needs a procedure to ensure that the monotonically decreasing curve contacts the horizontal axis at a definite finite time point t_{ p }. In other words, the cell population curve should be dynamically modulated by therapy so that its trajectory is actively guided, enabling the curve to contact the horizontal axis at time t_{ p }.
As in the case of a biochemical reactor, conditions may arise where the synthesis rate R following firstorder reaction, needs to be controlled and stopped at a definitive time [13]. To explore the possibility of inducing the reactor's processing to stop (i.e., R = 0), one may omit the input feed of reactants which causes the reaction to slow, but not become zero even after any length of time leaving a small asymptotic level of reaction rate R persisting at a time t, since firstorder kinetics indicates R′ + κ R = 0.
The asymptoticity in both the above cases is due to the fact that both the dynamic equations are similar and firstorder, viz. y′ + k y = 0. Here the yaxis denotes the system parameter, as reaction rate or altitude, while the xaxis denotes time. The equation implies that at infinite time, the curve meets the yaxis, if the coordinates at this meeting point is t^{†} and y^{†}, then y^{†} = 0, as t^{†} = ∞. Nevertheless, in the exponential systems above, we elucidate that the systems dynamics can be halted in a finite time (i.e., y = 0, at specific time t_{ p }), if one implements the concept of an adjustable or tunable bias shift, denoted by y*. This enables the performance curve to approach a value y → y* (a proposed predetermined negative value), as t → ∞ [Figure 2B]. This ensures that the curve trajectory intersects the time axis, and has an exact value y = 0 robustly, at a definite time t_{ p } [point P, in Figure 2B]. We have earlier methodologically simulated and validated realistically the procedure of bias shift, while respecting the requisite bounds or constraints imposed on the system [14]; thereby the exponential curve was induced to become zero within specific time interval, the error in the validation process being within 2.5%.
Using bias shift as baseline regulating principle
The aim of control systems design is to enable the tumour cell population (state T) to track and approach the T* value i.e. T → T*, where T* has a chosen bias shift whose value depends on the desired duration of tumour elimination, t_{ p }. From Figure 2B, the curve T is the arc path ACPD.
To enable the T curve to approach towards the T* value (bias shifted line L in Figure 2B), we use the guided tracking principle of control systems analysis [15]. Actually, the tumour cell population T follows the curve from point A to P in Figure 2B, but when the T trajectory reaches P (i.e., when population T attains zero value), the curve becomes horizontal, as all the tumour cells have become extinct, and thereafter, the trajectory will follow the axisline PEF onwards, where T remains zero throughout as time elapses. The T population does not diminish further towards a negative value as that is a biologically impossible feat. Further, we can treat the benign state of the subject (without any tumour cells) as the target baseline condition, and hence the undesirable tumour cell population state T is considered as the tumour system deviation or error E_{ T } from the benign condition. We define the error (Figure 2B) as E_{ T } = (T – T*), and the error’s time derivative as E_{ T′ } = (T′ – T*′). Thus eq. (3) becomes the errorequation which is calculated as:
Eq. (5) is derived by substituting the expressions of E_{ T } and E_{ T′ } above, in the errorequation. In eq. (5), the term T*′, i.e. the temporal rate of change of T*, is zero, since value of T* is a prefixed constant (as decided by the experimenter). Thus, dropping the term T*′ in eq. (5) furnishes the required tumour cell reduction rate for tumour elimination:
Eq. (6) is the condition for complete tumour regression, and hence should be followed by the tumour cell compartment in Figure 1. Thereafter, proceeding upwards through the successive compartments, such as administration of chemotherapy then the immunotherapy, we can obtain the dosetime profiles of temozolomide, interleukin2 and tumourinfiltrating lymphocyte injections, for regression of the tumour within the time duration t_{ d }. Thereafter, the tumour shall not recur as there are no surviving tumour cells. As a precaution, we continue the therapies for an extra sufficiently long period before fully stopping.
We have earlier used the inverse solution approach, but without bias shifting, for two treatment scenarios: (i) controlling chemotherapy infusion (imatinib) in myeloid leukemia [16], and (ii) regulating therapeutic infusion for treatment of ionic metabolic or hypocalceamic imbalance (control target error < 5%) [17]. In these cases the tumour cell population regressed and the blood ionic level tended to approach the desired value with high stability and asymptotically, but complete tumour elimination could not be obtained at a definitive time. In the present paper, we have remedied this problem by using adjustable bias shifting.
Computational model of multimodal therapy
We elucidate the quantitative approach for tumour dynamics under combinational chemotherapy and immunotherapy (Figure 3), based on coupled ordinary differential equations (ODE) between various cells and therapeutic agents. The approach is a broad one, and can be generalized to common situations. In the current instance, we take the case involving tumour cells, innate immunity cells (Natural killer cells), acquired immunity cells (Cytotoxic Tcells), and normal tissue protective cells (Circulating lymphocytes), along with antitumour therapies, including chemotherapy in the form of TMZ, and immunotherapy, namely externally administered lymphocytes (tumourinfiltrating lymphocytes) and lymphocyteactivating agents (as IL2). The equations recollect the experimentallybased immunological reactionrate framework demarcated by de Pillis et al.[18], Kuznetsov et al.[9] and Kirschner et al.[10], which has also been empirically validated [19]. The formulation also utilizes findings from animal studies and human clinical trials and uses models employed in cellular population biology, incorporating reaction kinetics such as MichaelisMenten type of interaction and saturation processes as logistic functions.
All the cellcell interactions, rate equations with temporal dynamics of the proposed model are graphically summarized in Figure 4 [eq. (S.1) to eq. (S.7) therein], and these are elucidated in Additional file 1: Supporting Methods. In the formulas, the primed symbols denote the temporal derivatives of the above entities (populations of the four types of cells, or concentration of three therapeutic agents). The numerical values of the constant parameters are obtained from experimental or clinical studies; values are given in Table 1.
Numerical values of the variables from the patient’s clinical data
Regarding the variables in the left side of eqs. (S.1)(S.7), N is the NK cell population, C is circulating lymphocyte population, L is the cytotoxic T cell population in blood, M is the chemotherapeutic drug concentration (in serum), I is the immunotherapeutic drug concentration (in serum), and T is the tumour cell population in the tumour. One may note that all these variables have definitive numerical values, which vary as time elapses and therapy continues. These values can be calculated by solving the equations and using the numerical values of parameters in the left side of eqs. (S.1)(S.6). Some of these parameters are related to the tumour cellular system:

(i)
The pharmacological/cell birth/death/interaction parameters: a, b, c, d, e, f, p, m, q, u, r _{ 1 } , r _{ 2 }

(ii)
Drug input/output parameters: α, β, γ, μ

(iii)
MichaelisMenten or saturation parameters: g, h, j, k, p _{ I } , g _{ I } , D, d, l, s;

(iv)
Temporal lysis parameters of different cells: k _{ T }, k _{ N }, k _{ L }, k _{ C }.
The other parameters of the equations are related directly to the therapeutic agents:

(a)
v _{M} (t), v _{I} (t), v _{L} (t) are the daily injected dosages (input rate per day) of chemotherapy (temozolomide), immunotherapy (interleukin2), and cytotherapy (TIL cells) respectively.

(b)
d, the saturation level of fractional tumour cell kill by the cytotoxic T cells,

(c)
l, the powerlaw exponent of fractional tumour cell kill by CD8+ T cells,

(d)
s, the steepness coefficient of the cytotoxic T cell – tumour cell interaction curve
Stationary condition of the system in tumour elimination
The equilibrium points or stationary conditions of the system can be found by setting the derivatives of the variables to zero, namely dT/dt = 0, dN/dt = 0, dL/dt = 0, dC/dt = 0, in eqs. (S.1)(S.4) respectively. Solving for the case t = 0, the equilibrium points are as follows:
This equation implies that there is a definitive tumourfree stationary state. At this equilibrium point there is no presence of any malignant cell, i.e., the cytotoxic Tcell population T_{E} = 0, L_{E} = 0, while there are specific levels of the immune cells, NK cells and circulating lymphocytes, that correlate with the healthy state (N_{E} = eα/βf, C_{E} = α/β). Substituting the values of the parameters α, β, e and f from Table 1 and taking the average blood volume of 6 liters for an adult, we arrive at the natural killer cell equilibrium population = 430,000 cells for the whole subject and circulating lymphocytes population = 72 billion for the same. It is this equilibrium state that the system tends to settle to, in the occasion of successful therapy.
System design formulation of multimodal therapy
In order to describe a system in terms of a control systems schema, we need to define the internal parameters (state variables characterizing the condition of the tumour system) and the external parameters (control parameters that can alter the tumour system). We formulate the ordinary differential equations (Figure 4), in terms of the control theory, whereby

(i)
System state variables are the parameters indicating the internal biological condition of hosttumour interaction in the patient. The system state variables in this case are (1) Tumour cells (T), (2) Natural killer cells (N), (3) Circulating lymphocytes (C), (4) Cytotoxic TLymphocytes (L), as well as (5) blood concentrations of Interleukin (I) and of (6) temozolomide chemotherapy (M) [Figure 1]. The state variables are denoted as functions X _{ n }, where n = 1 to 6.

(ii)
Control variables of system include the parameter values the external administrations of the three therapeutic agents, viz. the dose input rates of the injected chemotherapy (temozolomide) and immunomodulants (interleukin2 and tumourinfiltrating lymphocytes), which are represented by v _{M} (t), v _{I} (t), and v _{L} (t) respectively. These parameters vary with time and are measured in bodyweight normalized rate units, e.g. in mg. of a drug per kg of body weight (or, per square metre of body surface) of the patient, given per day. These control variables are denoted as functions U _{ m }, where m = 1 to 3.

(iii)
Constraint conditions of system indicate the quantitative requirements or thresholds that cannot be crossed application of the therapy. All the state variables need to be maintained within minimum and maximum values given in Table 2. For instance, the minimum and maximum values of the different cellular parameters should be within the physiological range that is necessary for homeostatic maintenance of the internal environment, or milieu interne of the patient. Too high a value of an entity can be toxic to the system, while, in some cases, there is a minimum value required so that the host can have proper immune surveillance to ward off foreign cells or microorganisms.
Using control systems practice, to construct the rate equation of any one particular state variable X_{ p } (a biological parameter), we express its time derivative or rate parameter X_{ p }′, as a function of (i) the values of the various biological parameters that determine the system, namely X_{ n }, and (ii) the effect of the proximal causal control variables (drug factors, U_{ m }) on the biological variables (X_{ n }).
We shall now explain the construction of the control design, in terms of objective, operation and toxicityminimization for each compartment.
Level 1: tumour cell compartment
Objective
The aim is to find the desired values, M* and L* of this compartment’s input parameters (Temozolomide blood level and Cytotoxic Tcell population), which will drive the compartment’s output (tumour cell population curve T) towards the line T*, until the T value becomes zero at point P in time t_{ P } (Figure 2B).
Operation of the compartment
For the tumour cell compartment in Figure 1, we construct the rate equation of the output state variable, tumour cell population T, by expressing its time derivative T', in terms of

(i)
the function f _{ T } (X _{ n }) which denotes the values of the various biological parameters that determine the tumour cell population system, namely the state parameters X _{ n. } This explains the tumour cell alteration function due to two intrinsic biological (nonpharmacological) processes: (a) tumour cell growth term, due to logistic or saturable growth, (b) tumour cell elimination term due to natural killer cells.

(ii)
the functions g _{T 1} and g _{T 2} that imply the biological effect of the compartment’s input therapy variables that show a dose–response saturation behavior U _{ r } , which is produced by the two antitumour entities, cytotoxic Tcell population L, and temozolomide concentration in blood, M.
Thus,
where U_{ M } and U_{ L } indicate the antitumour behavior or therapeutic efficiency of Temozolomide and cytotoxic Tlymphocytes respectively, while b_{ T } denotes the combined antitumour effect of these two therapeutic entities on tumour growth rate T′ (Table 3). This antitumour behavior is quantitatively a saturation type of dose–response function (like MichaelisMenten function, Hill function or similar activity function), which is a ubiquitous property of pharmacological responses. Thus Eq. (7) which is equivalent to eq. (S.1) of Figure 4, can be written in terms of MichaleisMenten kinetics as,
Comparing the right side of eq. (7) and (8),
and
Now, for enforcing the tumour elimination using the tracking approach mentioned earlier, we use the previously described tumour growth rate as T′ = −κ_{ T } (T – T*). [One may take care to distinguish this kappa parameter, κ_{ T, } the tumour decay rate (due to multimodal therapy), from the cytolysis rate parameter k_{T} of eq. (S.1) due to chemotherapy only]. Putting this expression of T′ = −κ_{ T } (T – T*) in eq. (7), we obtain the condition that needs to be followed to induce complete tumour regression:
Rewriting eq. (11), we obtain the terms containing the pharmacological effect functions g_{T 1} and g_{T 2}:
where,
Eq. (12) furnishes the relationship of the blood concentration of the antitumour agents (U_{ M } and U_{ L }), which, if implemented, will ensure that the tumour cells undergo complete extinction.
Transposing the eq. (13),
Eq. (13) implies that, for proper therapy, the combined therapyinduced antitumour effect term, b_{ T } , should take care of two aspects: (i) the condition required for enforcing tumour cell population elimination, i.e. the expression κ_{ T } (T – T*), and (ii) the tumour cell growth due to biological processes, i.e. f_{ T } (Table 3). At each time point, as tumour cell population changes, the cytotoxic therapy drive term b_{ T } alters. Eq. (13) indicates that the relaxation decay term b_{ T } has negative value if tumour is to regress (or b_{ T } has value zero, if the tumour becomes arrested at a definitive volume and thereby becomes stable), thus here b_{ T } ≤ 0. From eq. (12):
Eq. (12) furnishes the relationship of the blood concentration of the antitumour agents (U_{ M } and U_{ L }), which, if implemented, will ensure the tumour cells undergo complete extinction.
Minimization of toxicity of antitumour therapy
However the blood concentrations U_{ M } and U_{ L } also needs to be regulated in such a way so that the combined toxicity effect on the patient is minimal. Hence we now need to minimize the normal tissue damage cost function due to the agents. Using the quadratic cost concept (see earlier “Background” section), we have the cost function for the tumour cell compartment as:
Here r_{T 1} and r_{T 2} are the weighting factors of the two therapeutic entities, temozolomide and cytotoxic Tcells. The cost function J_{ T } needs to be minimized, and at the same time eq. (12) should be obeyed as a constraint. The minimization of the cost function, J_{ T } = ½ [r_{T1} U_{M} ^{2} + r_{T2} U_{L} ^{2}], is evidently a constrained optimization problem, and the optimality can be customarily solved by the method of Lagrange multiplier, λ. Thereby the augmented cost function J can be written as:
where the right sided expression within the second brackets {….} in eq. (16) incorporates the left side of the constraint equation [eq. 14(A)], as required by the Lagrange’s method. To pursue the minimization of the augmented cost function, we differentiate J with respect to the two variables U_{M} and U_{L}. Thereby, we find out the minimization conditions, i.e., ∂J/∂U_{ M } = 0, and ∂J/∂U_{ L } = 0, whence we get the expressions for the two variables:
Temozolomide efficiency term,
Cytotoxic Tlymphocyte efficiency term,
where,
Further, from eq. (10):
Desired input values of the therapy levels
Solving the two equations of the last line, we obtain the desired values, M* (blood concentration of temozolomide) and L* (cytotoxic Tcell population):
In the last two equations, for the terms U_{M} and U_{L}, we substitute their values from eq. (17) and (18) respectively. Thence, we arrive at the desired values of the temozolomide blood level and cytotoxic Tcell population, which, if implemented, will regress the tumour fully:
Desired temozolomide blood level:
Desired cytotoxic Tcell population:
Level 2: cytotoxic Tcell compartment
Objective
Here the goal is to find the desired values, I* and v_{ L }*, of this compartment’s input parameters (interleukin2 blood level and tumourinfiltrating lymphocyte injection doserate), which would drive the compartments output, namely the cytotoxic Tcell population, L to its desired value L* mentioned in the last paragraph. This enforced driving needs to be faster than the earlier compartment (tumour cell compartment) and requires to be done till time t_{r}, point P, when all the tumour cells have become eliminated (Figure 2B).
Operation of this compartment
We use the methodology of guided tracking principle of control systems analysis, which we have already used for the tumour cell population compartment [eq. (5)]. Likewise, we enforce the situation of the cytotoxic Tcell population L obeying the condition E_{ L }′ + κ_{ L }E_{ L } = 0, where κ_{ L } is a positive parameter and the deviate E_{ L } = (L – L*). Using similar logic of the reduction of the deviate or error [eq. (6)], we obtain:
We recollect the characteristic rate formulation of the tumour cell compartment, T′ = f_{ T } (X_{n}) + b_{ T } [i.e. eq. (7)]. We now desire to formulate the characteristic rate formulation of the cytotoxic Tcell compartment in terms of L′. Thus, we express the temporal dynamics equation of L′ [Eq. S.3 of Figure 4)] as:
where v_{ L } is the tumourinfiltrating lymphocyte injection doserate, U_{I} is the therapeutic efficiency relationship of interleukin2, and f(X_{n}) is a system function expressing the Tcell immunomodulation, viz. the biological and pharmacological variables, acting on this cytotoxic Tcell population. Compare eq. (25) with eq. (S.3) that states:
By this comparison, one notes that U_{ I }, the term dealing with interleukin level I, can be identified as:
and one can also express
Eliminating L′ from eq. (24) and eq. (25), we note that
Where
Note here that f_{ L }(X) is the performance function of the cytotoxic Tcell compartment, a part of the the rightside of the eq. (26) except its last two terms, which are the therapyinput terms, dependent on the interleukin and tumourinfiltrating lymphocyte injected, the two terms being denoted as b_{ L }, In eq. (26), the last term v_{L} is the Tumourinfiltrating lymphocyte injection doserate, while the secondlast term [p_{ I }LI/(g_{ I }+I)], now denoted as U_{I}, is the therapeutic efficiency relationship of interleukin2. Evidently from eq. (26), the term b_{ L } signifies the total cytotoxic Tcell activation effect by the two immunotherapeutic inputs: the interleukin efficiency term U_{ I } and the tumourinfiltrating lymphocyte administration term v_{L} (Table 3). Indeed, the said equation indicates that the relaxation decay term b_{ L } has positive value if the tumour regresses (or zero, if the tumour is arrested and thereby stable), i.e. b_{ T } ≥ 0. Actually, eq. (29) furnishes the relationship of the characteristics of the antitumour agents, interleukin and tumour infiltrating lymphocytes (U_{I} and v_{L}), which if implemented, will ensure the full elimination of tumour.
Minimization of toxicity of antitumour therapy
In this compartment, the two antitumour agents are nterleukin2 and Tumourinfiltrating lymphocytes. As explained for the earlier compartment, we need to minimize the toxicity cost function for the compartment where the two therapeutic input functions arriving at the compartment are Interleukin2 (i.e., its efficiency function U_{I}) and tumourinfiltrating lymphocyte (its injection doserate, v_{ L }). The cost function for this compartment is:
where r_{L 1} and r_{L 2} are the sensitivity weights due to the two aforesaid therapeutic moieties in this compartment. For minimization, the constraint requirement [U_{ I } + v_{ L } = b_{ L }] from [eq. (29)] should be obeyed. Solving using the Lagrange multiplier method, we arrive at:.
Interleukin2 efficiency term,
Tumourinfiltrating lymphocyte injection doserate,
where
Desired input values of the therapy level
By transposing U_{ I } = {(p_{ I }L_{ I })/(g_{ I } + I)} [i.e., eq. (27)], we obtain the desired value of I* (blood concentration of Interleukin2):
Substituting U_{I} from eq. (32) to eq. (35), we get the desired value of tumourinfiltrating lymphocyte injection doserate which will eliminate the tumour:
Desired interleukin2 blood concentration,
Level 3: temozolomide injection compartment
Objective
Here the goal is to find the desired value, v_{ M }*, of this compartment’s input parameter (temozolomide injection doserate), which would drive the compartment’s output, namely the patient’s blood level of temozolomide M to its desired value M* as given in eq. (22). This driving needs to be done faster than the preceding compartment (Tcell compartment), and is to be done by time t_{ p }, point P (Figure 1).
Operation of this compartment
We need to implement the condition that E_{ M }′ + κ_{ M }E_{ M } = 0, where the deviate E_{ M } = (M – M*). Using similar reasoning as the earlier compartment, we get M′ = −κ_{ M } (M – M*). Eliminating M′ from this equation and from eq. (S.5), we have:
Level 4: interleukin2 injection compartment
The condition necessary to be implemented is that E_{ I }′ + κ_{ I }E_{ I } = 0 where the deviate E_{ I } = (I – I*) Thus, I′ = −κ_{ I }(I – I*). By eliminating I′ using eq. (S.6), we note:
Level 5: tumourinfiltrating lymphocyte injection compartment
For this compartment, we can write from eq. (33) that:
Determining tumour regression rate constant, bias shift and therapeutic weights
Delineating the rate parameters κ_{T}, κ_{M}, κ_{L}, κ_{I} and bias T*
These are rate constants of the tumour cell compartment, temozolomide compartment, cytotoxic Tcell compartment and interleukin compartment, respectively, as regression occurs under the action of multimodal therapy. These are calculated from the desired rate of tumour regression, expressed as settling time t_{s} of the regression process (the time duration in which 90% of tumour has regressed), and is taken to be around 1–2 months. The tumour regression rate constant κ_{ T } = 4/t_{s}; so if t_{s} is 60 days, κ_{T} = 0.067 per day. On the other hand, the dynamics of the successively preceding compartments (e.g. chemotherapy and cytotoxic Tcell compartments) need to be faster, as they causally influence the tumour cell compartment, and thus need to change more rapidly if they are to have a controlling influence on the tumour cell compartment (Figure 1). Thus, the time constants of the modular stages will be lesser, and hence the process rate constant will be higher. So, we can choose the rate constants κ_{ L } and κ_{ M } of cytotoxic Tcell compartment and temozolomide injection compartment respectively, such that they exceed κ_{ T }. Similarly, the rate constant κ_{ I } of interleukin compartment (that causally acts on the Tcell compartment), is chosen to be higher that κ_{ L }.
Thus,
As explained earlier, in control analysis [15], a thumb rule used is: rate constants between two successive stages in a cascade control system are in the ratio of 1: 3. Thus, if we have the value of rate constant of the tumour cell population, κ_{ T } and then we obtain the value of the other constants in terms of κ_{ T } :
In the above example, since, κ_{ T } = 0.067/day, we get κ_{ L } = 0.201/day, κ_{ M } = 0.201/day, and κ_{ I } = 0.603/day. It may be noted that the tumour elimination time t_{ p } (Figure 2) when 100% tumour has regressed, is longer than t_{s}, and can be selected as 10 days more, i.e. t_{ p } = 70 days. If the initial pretherapy tumour cell population is estimated as T_{ 0 } , then substituting the values of T_{ 0 }, κ_{ T } and t_{ p } in the rightsided expression in eq. (4), provides the value of bias shift T*.
Selection of therapeutic weights r_{T1} , r_{T2} and r_{L1} , r_{L2}
These adaptable parameters, r_{ T1 }, r_{ T2 }, r_{ L1 }, r_{ L2 } are needed for minimizing the toxicity cost of the therapy. It is these parameters that give a control to the investigator for maneuvering the tumour elimination process, under adjustable dosing of the drugs. Initially, the values of the tuning parameters r_{ T1 }, r_{ T2 }, r_{ L1 }, r_{ L2 } which appear in derivation, are chosen by specific quantitative conditions (see Additional file 2: Supporting Analysis). To recapitulate, r_{ T1 } and r_{ T2 } are respectively the toxicity cost weighting factors of temozolomide and cytotoxicTcell, acting on the tumour cell population compartment, producing the cost J_{T} which is to be optimized.
For such optimization problems, it is known that the important characteristic to be considered is the ratio r_{ T1 }: r_{ T2 }[15]. Thus, we can take the parameter r_{ T1 } to have a normalized value of unity (i.e. r_{ T1 } = 1), thereby the task is to suitably choose or optimize the value of the other tuning parameter r_{ T2 }. Correspondingly, r_{ L1 } and r_{ L2 } are respectively the toxicity cost weighting factors of interleukin2 and tumourinfiltrating lymphocytes, which act on the cytotoxic Tcell compartment and produces the toxicity cost J_{L}. Similarly, this cost can be minimized by normalizing r_{ L1 } = 1, and then optimizing r_{ L1 }. All the tuning parameter values must be greater than or equal to zero. Using the upper and lower bounds that needs to be followed by the cellular and pharmacological variables (Table 2), we can calculate the numerical values of the therapeutic weight parameters (Additional file 2: item B.3 and Table 1 there).
Efficiency of different combinations of therapeutic agents
During the analysis and modeling procedure, the control variables, i.e. parameters of the therapeutic agents, should have positive values and not attain negative nor imaginary values, even though the system dynamics [eq. (S.1S.7)] has fractional exponent at some places [eq. (S.7)]. The situation is amended in case if any of the computed control variables, or of the three injected drug dosage rates calculated (v_{M}*, v_{I}* or v_{L}*), becomes imaginary or negative. If this happens, the particular drug injection rate (say, for drugA) is stopped and kept suspended (its dosage is henceforth taken to be zero), the simulation is continued with the presence of the other two drugs. As long as the dosage of the drugA calculated remains negative or imaginary, drugA is not injected, and the drug A’s value for the subsequent simulation steps becomes zero, while the other two drugs, whose calculated dosages continue to be positive numbers, are continued to be administered. Only when the computed dose of the drugA turns positive, then the injection of drugA is resumed according to that dose. To wit, the following alternative situations occur:

(i)
The dosage parameters of all three drugs remain positive throughout the time duration: If so, all are given.

(ii)
One drug becomes negative or imaginary for a particular time interval: this drug is stopped during that interval, and the other two drugs are continued as per the modeling procedure with the first omitted

(iii)
Two drugs become negative or imaginary: here the remaining drug with a positive drug dosage rate is administered.

(iv)
All the three drugs violate the positivity condition: then no drug is given, until one or more drug rates become positive at a subsequent time, upon which the administration of the drug/s is resumed.
If case (ii) is applicable, then there will be three different combinations for two drugs at a time by permutation: (a) temozolomide and interleukin, or (ii) temozolomide and tumour infiltrating lymphocytes, or (iii) interleukin and tumour infiltrating lymphocytes. If case (iii) is relevant, then the three possibilities are: temozolomide or interleukin or tumour infiltrating lymphocytes. The mathematical formulation for cases (ii) and (iii) are derived later (Additional file 2: item B.4), using the overall approach of case (i) presented above. Further, once the tumour cell count becomes zero (point P, Figure 2B) by any of the above approaches, one stops the drug/s thereafter. The overall methodology of our approach is shown in Figure 5.
Results and discussion
Numerical simulations for tumour elimination
We computationally model the tumour system dynamics using custommade coding and solve in the MATLAB platform [eq. (7)(14A)], and thereby obtain the temporal profile of the varying doserate of the therapeutic agents, by solving eq. (37)(39) which respectively furnish the temporal profile of temozolomide, interleukin2 and tumourinfiltrating lymphocytes, which would eliminate the tumour cell population. In Figure 6 we develop the detailed procedural steps into a readilyusable algorithm toolkit. Our concern is the general problem of neurooncology, especially neuroectoderm originating tumours such as glioma, melanoma or neuroblastoma, which are comparable tumours, biologically and pharmacologically, and respond to similar therapeutic interventions, as temozolomide, interleukin2 and lymphocyte immunotherapy by Tcells [1, 5, 20, 32, 33]. Now we apply the proposed approach to regressing both slowly and rapidly growing tumours, astrocytoma grade II and astocytoma grade III respectively, aiming at a timeframe of about 2 months. We take the values of parameters of the equations from experimental clinical studies available or from general oncological investigations (Table 1). The values of the parameters can be applied to tumours in general, for example to tumours of neuroepithelial origin, as melanoma or neuroma, or sarcoma [9, 19, 24]. Actually, using values of such parameters from clinical or preclinical tumour model systems, one can solve differential equations of tumour response under drugs or immune cells, and the mathematical predictions are closely confirmed and validated by cytological measurements on human subjects during therapy [19, 34, 35].
Lowgrade tumour
We consider the realistic case of lowgrade tumour as ologidendroglioma or glioma grade II tumour with initial conditions of starting malignant cell population (T_{ 0 }) of 2 × 10^{7}, natural killer cell population of 10^{5}, cytotoxic Tcell population of 5 × 10^{4} and circulating lymphocytes of 10^{9} cells (see endnote ^{b} before the References). From this T_{ 0 } value, we get a bias shift T* = 1.8986 × 10^{5} cells [rightside of eq. (4)] (see subsection on “Determining tumour regression rate constants” above). All the other constants used in the model are given in Table 1, with the tumour cell growth rate, a = 0.301 per day, and the deceleration rate of logistic tumour growth b = 1.01x10^{8}. We start the simulation procedure increasing the temporal step by one minute each time (Figure 7), and find that the tumour cell population actually become zero for extinction at 30.1 days (Figure 7A). One may notice that in this case the tumour cells are eliminated before the desired time duration, t_{ p }. This duration is different from the initially planned tumour elimination time t_{ p }, since, during simulation there may be time durations when a drug is stopped [this happens if its dose is calculated to have a negative or imaginary value [Figure 6, item (M)]], but the other drugs are increased with sufficient intensity to ensure tumour cell lysis at a substantial level. These increased therapeutic dosages can have nonlinear cellkilling effects and thus produce a shorter tumour elimination time than expected. It may be noted that we perform the computation for a time span of 20 times the extinction period (i.e. span of 800 days), a sufficiently long time to check permanency of regression. The regression appears to be lasting without any reappearance of tumour cells, even though the therapy had been stopped much earlier (Figure 7A). It may be noted that for our plots or records, we take the computed cellular population values to have the next integer numeral, e.g. if N, C, T or T* is calculated out to be 1230.23 cells, the value is taken to be 1231 cells.
We also plot the blood concentrations of the cytotoxic Tcells, chemotherapeutic agent temozolomide and interleukin2, required to eliminate the tumour (Figure 7B7D). Also displayed are plots of the population of circulating lymphocytes and NK cells to check that there is no significant toxicity as sideeffects of the therapy (Figure 7E7F). We also observe that these values are also well within the corresponding bounds of the human host system (Table 2). Finally Figure 7G7I show the injected dose rates, as required, daywise, for each of the three therapeutic agents, that enables full tumour elimination. The circulating lymphocytes and NK cells takes around 500 days to reach the steady state values of 76 billion and 510,000 cells respectively. Note that these values closely corroborate with the stationary points solved theoretically earlier, which are 72 billion and 430,000 cells respectively [see eq. (6A) and subheading “Stationary condition of the system” there]. We find that the extreme values of these variables in the graphs are well within the bounds of the human host system (Table 2). For instance, the maximal doserate of temozolomide chemotherapy and of tumourinfiltrating lymphocytes are less than 10% and less than 1% respectively, of the upper bound of these agents in Table 2. Further, the cumulative dose of the interleukin (from the doserate graph calculated as area under the curve in Figure 7D), is below 1% of the interleukin upper bound in Table 2.
Highgrade tumour
For the faster growing tumour such as astrocytoma or glioma grade III, we take tumour cell growth rate, a = 0.431, which is about 50% higher than the lowgrade glial tumour discussed above. An initial malignant cell population (T_{ 0 }) of 2 × 10^{7}, natural killer cell population of 10^{5}, cytotoxic Tcell population of 5 × 10^{4} and circulating lymphocyte of 10^{9} cells were considered to simulate the system. The same tumour settling time and tumour elimination timeline is given, and a bias of T* = 1.8986 × 10^{5} cells is taken. The other parameters for the tumour system are given in Table 1. We proceed with the numerical simulation and observe that the tumour cells become extinct at 40.11 days (Figure 8A). Since the tumour cell proliferation rate, a value is higher here than the earlier tumour, the current tumour needs more treatment duration than the former. Likewise, we perform the computation for a time span of 800 days, to ensure the permanency of regression far across time. The values of the biological and pharmacological entities are displayed in Figure 8B8I. The circulating lymphocytes and natural killer cells stabilize around 0.33 million and 62 billion cells respectively (Figure 8E8F). These correspond to the stationary state populations predicted mathematically earlier. As in the earlier case of slowgrowing tumour, one can also note the same for the rapidlygrowing tumour:

(i)
the dosages of the therapeutic agents drops to zero from the tumour extinction time onwards

(ii)
none of the cell populations or dosing of therapeutic agents crosses the bounds (Table 2), and

(iii)
the dosages of the three agents are approximately the same fraction of the upper bound (10%, 1% and 1%).
We also discern that the therapy period is considerably longer in highgrade tumour, basically since the tumour growth rate is about 50% more than the slowgrowing tumour (value of parameter α, in Table 1). Further, since these therapeutic dosages are a much smaller fraction of the upper bounds, the dosage profile can be well tolerated by patients. Actually, the overall patterns of the graphs (Figures 7 and 8) are comparable, across both the high and low grade tumours.
Robustness of tumour elimination procedure
The values of the parameters used are taken from Table 1. However these biological system parameters can vary across individuals, depending on the physiological and constitutional of the patient. Hence we need to check that the methodology is applicable even though the parameters may vary along a range. In order to test the versatility and robustness of our approach, we varied the initial immunological status (values of N, C, L cells) and tumour system parameters (values of tumour growth rate, α and carrying capacity, β). We performed 500 random simulations to check for complete tumour elimination, in both tumour grades. We found that tumour extinction for oligodendroglioma occurred in 100% cases if the variation was 0%, and in 98% of the cases as the coefficient of variation was increased to 10% (Table 4). For astrocytoma, the results are comparable. It has been known that the physiological parameters are generally kept constant homeostatically by organisms, with a 10% variation around the mean level [31]. Thus we can ascertain that the proposed approach may be able to induce tumour elimination in the large majority of cases.
Unitary pattern in tumour regression process
If one compares the corresponding graphs in Figures 7 and 8, there are evident similarities in the temporal profiles of the therapeutic agents needed in both highgrade and lowgrade tumours to enable elimination of malignancy. The common patterns valid across both tumours are elucidated below.
Profile A: terminal therapy pulse and cytotoxic lymphocyte persistence
One notes that (i) no injections of any of the three therapeutic agents are required any time after the time point of extinction of tumour, which does not recur later (Figure 7G7I and Figure 8G8I); (ii) the blood levels of the three therapeutic agents and the populations of the natural killer cells and circulating lymphocyte stay within the requisite limits, assuring that there would be no significant toxicity to the patient (Figure 7B7F and Figure 8B8F). One can make two pertinent observations from the simulation results. Firstly, for one to make the exponential decay curve of tumour cell population hit the T = 0 baseline (i.e. the xaxis of Figure 2B) at a definitive time point t_{ p } , one needs to inject a terminal pulse of each of three therapeutic dosages before the end of treatment protocol in both the lowgrade and highgrade tumours, as shown by arrows in the six graphs in Figure 7G7I and Figure 8G8I. The conjoint effect of the pulses of all three agents ensure that in the last stage, all the tumour cells in the exponential extremity of cell population decay curve, do become eliminated. Secondly, the persistence of the cytotoxic T lymphocytes in blood for a over a week after tumour has been eliminated (Figures 7B and 8B), can have a beneficial effect, such that it can act as a vigilant antitumour measure against recurrence, by acting for an appreciable time after tumour extinction. After that duration elapses, this lymphocyte population tapers off.
Profile B: common temporal paradigm of therapeutic agents needed
An insight into the mechanism of tumour regression may be obtained by investigating the commonalities in the pattern of the temporal variation in the blood levels of the therapeutic agents and population of the protective cells, which seem to be common across both highgrade and lowgrade tumours. From Figure 7B7D and Figure 8B8D, we discern the following temporal patterns of the aforesaid entities, the pattern being similar for both rapidly and slowly growing tumours:
(i) Temozolomide concentration (chemotherapy): monophasic activation pattern
The curve of M is unimodal for both lowgrade and highgrade tumours, steadily rising initially, and later, gradually tapering off (Figures 7C and 8C). This chemotherapy concentration curve follows the usual pharmacokinetic curve of a customary drug administration: a unimodal hump, with gradually rising drug concentration to point D, and, later gradual concentration reduction (Figure 9A shows overall pattern).
(ii) Cytotoxic Tlymphocyte (CTL) concentration (cytotherapy): biphasic activation pattern
This graph has bimodal peaks for both the highgrade and lowgrade tumour (Figures 7B and 8B, arrows; Figure 9B). We observe that for tumour elimination, there is a need to have a peak of CTL, both, around the initial and the final phases of the therapy. The initial peak concentration in required to forcefully target and guide the trajectory of the cancer cell population curve (T) towards the zero baseline (Figures 7A and 8A). The final peak concentration of CTL, occurs before the time of tumour elimination, and is necessary to sufficiently eliminate the tumour cells and depress their population trajectory so that the same hits the zero baseline at the definitive selected time point P (Figure 1B). Depending on the tumour system dynamics, the height or activationlevel of the second peak may be lower (Figure 7B), or higher (Figure 8B), than the height or activationlevel of first peak. In case of Cytotoxic Tcell concentration, the first peak is due to therapeutic agents of tumourinfiltrating lymphocytes and interleukin, whose dosing starts as an impulse stimulus from the initial time. The impulse of these agents also enhance the generation of CTL population. The second peak in cytotoxic Tcells (point B in Figure 8B) occurs due to the later dip or decay of chemotherapy concentration along the hump of the M curve (point K in Figure 8C). The chemotherapeutic agent is toxic to and diminishes all the cellular constituents, including CTLs, and, hence, a decrease of chemotherapy induces a rise of CTL then. This same pattern of primary and secondary induction of Tcell population also occurs in the other tumour (Figure 7B7C).
(iii) Interleukin2 concentration (immunotherapy): stationary activation pattern
This agent is needed at a substantially high level (at the upper bound, just below toxicity level), so as to induce a higher level of immune activation that would enable complete elimination of glioma cells (Figures 7D, 8D and 9C). Indeed it is well known from clinical experience in neurooncological immunotherapy that interleukin2 administered at significantly augmented dose, induces longlasting immunomodulation to act against those malignant cells that bypass usual therapeutic intervention [36]. Actually, the graph for the high and lowgrade tumours has rapidly rising high amplitude, as the interleukin2 level is truncated and kept stationary within toxicity limit.
(iv) Circulating lymphocyte population: saturating activation pattern
This population initially increases and then plateaus in both tumours (Figures 7F, 8F and 9D), to approach the saturation level of the longterm steady state as mentioned earlier.
(v) Natural killer cell population: saturating activation pattern
Similar to the Circulating lymphocytes, the NK cells show saturation behaviour for both high and low grade tumours (Figures 7E, 8E and 9D).
Experimental, biological and clinical corroboration
To justify, the modus operandi of our approach is based on the following aspects: firstly, the tumour elimination process is a manifestation of exponential diminution with bias shifting; secondly, the procedure works on fundamentally two complimentary modes of tumour cell reduction:

(a)
Reducing the tumour cell proliferation, as by chemotherapy or chemical alkylation (DNA damage);

(b)
Increasing the tumour cell lysis by antitumour lymphocytes, which are activated by cytokine modulation (for instance tumour infiltrating lymphocytes, interleukin2 etc.).
Further, the model developed shows that to induce tumour regression, the three antitumour entities should have three distinct temporal profiles: (1) biphasic intensity for lymphocyte activation, (2) monophasic intensity for activation of chemomodulative DNA damage (chemical alkylation), (3) stationary intensity of cytokine activation (interleukin2). The experimental studies in animals [37–39] and clinical situation [40, 41] corroborate these findings. Thus our results show that complete elimination of tumour can be attained by (i) the fivepattern profile: activation of antitumour lymphocyte (bimodality), chemomodulative DNA damage (unimodality), interleukin (stationary), natural killer cell (saturation), and circulating lymphocyte (saturation), (ii) the two kinetic conditions: bias shift and exponential decrementing dynamics [the tumour trajectory formula of eq. (4)].
Translational applicability
There are two aspects where the approach could be improved. Firstly, one can increase the system robustness, which diminishes as the individual patientspecific fluctuations increases (see earlier subsection on Robustness). For realtime implementation, we can use the neuroadaptive controller, which can reliably follow the desired mathematical trajectory of the tumour cell population curve (Figure 2B), and can be well adapted to different values and fluctuations of biological parameters of different patients. We have used such a controller in devising imatinib chemotherapy dosing in chronic myeloid leukaemia [16], and 100% robustness was obtained even when the maximal tumour density varied from 150,000 to 400,000 cells/mm^{3}, indicating about 250% variation on the baseline level. Secondly, for proper monitoring of a patient, the declining tumour load T can be weekly or semiweekly estimated noninvasively, by MRI amideproton transfer imaging, which maps the cell proliferation intensity by amide mobility [42]. This emerging technology holds high potential, as this method is an efficient one to distinguish between various oncological conditions such as tumour recurrence, tumour haemorrhage or tumour necrosis.
In Figure 7A the tumour cell diminution trajectory is shown for the subject in question, the cell population being mathematically derived given the physiological parameters. It would be preferable if one could check this mathematically obtained curve T(t) with the actual reallife tumour cell population T (t) as the treatment process goes over months, the actual population T (t) being estimated radiologically at suitable discrete time intervals of τ (Figure 10). The aim is to have correctionmaking feedback so that the actual realistic curve T tracks and follows the desired computationallyoptimized curve T. In case if there is a discrepancy at a particular time point t_{ k }, between the actuallymeasured value and the mathematicallypredicted value (i.e. between T and T at time point t_{ k }), then the value at that moment of the actual number of tumour cells T_{ k } can be used as the input value of tumour cell population at time point t_{ k } in the simulation procedure for going to the next time point t_{ k+1 } as per the procedure in Figure 6. Thereby, the new dosages of the drugs can be calculated for reaching the same target, namely zero population of tumour cells. The result is that the actual trajectory of the diminution of the tumour cell population becomes a discrete approximation of the continuous mathematical curve, with the time interval τ (Figure 10). Nevertheless, the final target of zero tumour cell population is attained.
For monitoring while the therapy is in progress, one need to have at suitable time intervals, the overall measurements of the serum temozolomide and interleukin2 levels, as well as the populations of cytotoxic Tcells, natural killer cells and circulating lymphocyte in the blood (Figure 10). Only approximate values of these parameters or of the tumour cell population are required, as even large variations of these parameters can be accommodated by the neuroadaptive controller mentioned above. During practical implementation, it may not be possible to measure all the biological parameters, neither a parameter can be estimated at all the desired time points as there may be discontinuous or missed sessions. Here, one can use suitable quantitative procedures as the Kalman’s technique or constraint filter method which has been used in immunology to estimate the likely values of biological parameters at the discrete or missed sessions, given their measured values in other sessions [43, 44]. Thus, it would suffice to monitor the tumour parameters weekly or twiceweekly, to enable the therapy control system tracking and adapting function, as demarcated in Figure 10.
One may observe that during therapy (Figures 7E7F and 8E8F), the circulating lymphocytes and natural killer cells increase but there is no hazard, as the populations are within the upper physiological limits of Table 2. Moreover, for making the methodology more suited to personalized medicine so that we can preselect the most sensitive drugs beforehand, we can utilize the emerging method of tumourgraft technology whereby one can grow a therapeutically faithful model of the individual patient’s tumour biopsy tissue on experimental mice [45]. Tumour graft platforms can test different drug combinations, and preselect the most sensitive drug before starting the treatment, the accuracy of selecting effective drug molecule being 86%. Lymphocytic (CTL) immunomodulation can be of tactical utility, as it exhibits a range of unique behaviors [46], that chemotherapeutic drugs cannot, such as (i) the cells can migrate to the antigenbearing primary or secondary growths of tumour, even in hidden tissue depths, (ii) CTLs can continue to multiply automatically in response to immunogenic proteins of malignant cells, until all those cells become extinct.
Conclusion
We have formulated the distinctive dosetime relationship of chemotherapy and of immunotherapy (interleukin2 and cytotoxic lymphocyte), such that their orchestrated functioning, that incorporates a biphasic temporal profile for Tcell, ensures that all the tumour cells are eliminated within a desired time. Excess toxicity to the host is avoided, as the circulating lymphocyte and natural killer cells in blood are protected. The approach is patientspecific as the formulation depends on the tumour load, the levels of cytotoxic T lymphocytes, cytokine interleukin, natural killer cells and circulating lymphocytes. All these parameters vary with the individual patients, and hence the different therapeutic dosetime profile obtained for each specific patient will be optimally suited from her. The formulation put forward use of an innovative approach of bias shift, control systems analysis, performance cost minimization and inverse construction of drug input. The limitation of the proposed model is that we have not considered the effect of tumour cells resistant to the chemotherapy drug and also not taken care of angiogenesis process where the drug permeability to the tumour cells will hinder. A broadbased gamut of findings from animal experimentation and clinical investigations are shown to corroborate with the tumour extinction approach developed. An unanticipated but noteworthy finding is the importance of giving a terminal pulse of the therapeutic agents before the end of the therapy, so that all the tumour cells become extinct and there is no extra druginduced toxicity, the natural killer cells and circulating lymphocytes being within the physiological limits. This is in contrast to the generally prevailing view in clinical medicine, which advocates tapering off of the therapy in the later stages. This tapering may cause tumour recurrence in clinical praxis, as there is no intensive spike of the therapeutic agents to eradicate all the malignant cells.
To summarize, information from temporal dynamics of both the endogenous and exogenous tumour regression has been used to explore the mechanism and elucidation of integrative functioning of various therapeutic modalities, whose combined effect eliminates the malignant tumour as corroborated with experimental findings. The method proposed can be of wideranging application, and can be adapted for application to conditions where there is involvement of chemotherapy and/or immunotherapy. The procedure delineated is also applicable to other tumour systems, as it offers a principled approach to tumour containment and thus an incisive prospect for probing towards further biological and clinical situations.
Endnotes
^{a} Clarifications on upper/lower limits of biological and therapeutic parameters (Table 2 ).
Circulating lymphocyte population: The normal blood volume under active circulation is 3.54.5 litres. The higher limit of lymphocytes in individuals can be up to 100 × 10^{9} cells/litre [26]. Taking the upper bound of the blood volume, this translates to the higher bound of lymphocyte in a person to be 4.5 × 10^{11} cells. In contrast, the lowest lymphocyte count for patients for duration of 6–12 months across therapy can be 663–1160 cells/μl [27]. Taking the lesser value of the cell count and the lower level of blood volume, we get minimal bound of lymphocytes as 2.32 × 10^{9} cells, which the patient can tolerate up to 6 months.
Natural Killer Cells population: The upper limit of NK cell is 13% of lymphocyte population, with CD56/CD16 surface protein being the marker for these lymphocytes [28]. As the higher bound of lymphocyte population in the earlier paragraph is 4.5 × 10^{11} cells, we have the maximum value of NK cell population in the individual to be 5.85 × 10^{10} cells. On the other hand, the lower limit of NK cells (CD56/CD16 lymphocytes) is 0, occurring in people having natural killer cell deficiency condition [29], and a period of 3½ months have been noted for elapse of this condition, before considerable infection can set in [29]. Hence we also mention this time duration for the NK cells lower bound.
Tumourspecific Cytotoxic (CD8+) Tcell population: The upper limit of this cell population for a patient is 20150 cells/mm^{3}[25]. Using a blood volume of 4.5 litres, the total Tcell (CD8) population will be 6.05 × 10^{10} cells. Furthermore, tumourspecific cytotoxic Tcells, that are specifically active against a particular malignant lesion, has been known to come into play if the tumour is present, and to decay away if the tumour undergoes elimination [30]. Hence, one mentions the lower bound of these T cells to be 0.
Temozolomide chemotherapy infusion dosage rate: Maximum daily dosage [1] of temozolomide is 200 mg/m^{2}/day, i.e., 4.45 mg/kg body. wt. per day. One may choose not to give it, so the lower limit is 0.
Interleukin2 immunotherapy infusion dosage rate: Maximum infusion given is 7.2 × 10^{4} International Units (I.U.)/kg/day [1]. The lower bound can be set to nil, as above.
Tumourinfiltrating lymphocytes (TIL) immunotherapy cumulative dosage: The maximum cumulative dose for a patient over the whole duration of therapy is 13.7 × 10^{10} cells [25]. Likewise, the lower limit is zero.
^{b} Determination of initial number of the cells for numerical experimentation.
Tumour cells: The number of cells in a tumour which becomes clinically detectable is 10^{8} occupying a volume of 1 cm^{3}, out of these 10^{7} cells are malignant, while the rest are stromal cells [47]. For our quantitative experimentation, we take double the amount of tumour cells to have a safety factor of 2, thus giving 2 × 10^{8} malignant cells as our initial condition.
Circulating lymphocytes: We also note the range of normal values: leucocyte count = 4000 to 11000/mm^{3}, the fraction of lymphocytes are 1540%, and actively pulsating blood volume under circulation is 3.54.5 litres. To be cautious for ensuring tumour regression, we will consider the lower values in the range. By multiplying the requisite aforesaid values, we have the circulating lymphocyte population in the person as 2.1 × 10^{9} cells. Again to be on safer side, we take half of the population for our study (i.e. 1.05 × 10^{9} cells, say 10^{9} cells).
Natural killer cells: The fraction of NK cells is 113% of circulating lymphocytes [28], we take the lower value for calculation, and use the circulating lymphocyte value of 1 × 10^{9} cells given above. Further, there is available data on the antitumour effector factor of NK cells during study of endogenous tumour regression, namely 1:50 as target tumour cell: effector NK cell ratio, i.e. a value of 2% [48]. Multiplying these values, we get the effective population of tumourtargeted NK cells as 2 × 10^{5} cells. As earlier, we will consider half this value to be on the safe side, i.e. NK cell population = 10^{5} cells.
Cytotoxic Tcells: The fraction of cytotoxic CD8+ Tcells is 13% of circulating lymphocytes [49]. The cytotoxicity factor of activated tumourinfiltrating cytotoxic T cells (CTL) is 9.5% as regression process is underway [48]. One also knows that if the tumour is invasive and spreads, the fatigue factor come into force, which can cause the CTL efficiency (as estimated by cytokine production) to decrease to 12.2% of the level as compared to when there is no invasion of tumour [50]. Using these factors, we arrive at the effective population of cytotoxic T cells as 1.15 × 10^{5} cells, say 10^{5} cells. As per the lower end, we have taken half, i.e. 5 × 10^{4} cells as the CTL population.
References
 1.
Perry M: The Chemotherapy Source Book. 2001, Philadelphia: Lippincott, 3
 2.
Murphy G: Lawrence W, Lenhard R: American Cancer Society Textbook of Clinical Oncology. 1995, American Cancer Society: Atlanta
 3.
Newlands ES, Blackledge GRP, Slack JA, Rustin GJS, Smith DB, Stuart NSA, Quarterman CP, Hoffman R, Stevens MFG, Brampton MH, Gibson AC: Phase I trial of temozolomide. Br J Cancer. 1992, 65: 287291. 10.1038/bjc.1992.57.
 4.
Wheeler C, Asha D, Gentao L, et al: Clinical responsiveness of glioblastoma multiforme to chemotherapy after vaccination. Clin Cancer Res. 2004, 10: 53165326. 10.1158/10780432.CCR040497.
 5.
Fenstermaker R, Ciesielski M: Immunotherapeutic Strategies for Malignant Glioma. Cancer Control. 2004, 11 (3): 181191.
 6.
Mustafa M, Buchanan G, Winick N: Immune recovery in children with malignancy after cessation of chemotherapy. J Pediatr Hematol Oncol. 1998, 20 (5): 451457. 10.1097/0004342619980900000008.
 7.
Nanda S, Moore H, Lenhart S: Optimal control of treatment in a mathematical model of chronic myelogenous leukemia. Math Biosci. 2007, 210: 143156. 10.1016/j.mbs.2007.05.003.
 8.
Kogan Y, Forys U, Shukron O, et al: Cellular immunotherapy for high grade glioma: Mathematical analysis deriving efficacious infusion rates based on patient requirements. SIAM Jour Appl Math. 2010, 70: 19531976. 10.1137/08073740X.
 9.
Kuznetsov V, Makalkin I, Taylor M, Perelson A: Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bull Math Biol. 1994, 56: 295321. 10.1007/BF02460644.
 10.
Kirschner D, Panetta J: Modeling immunotherapy of the tumor–immune interaction. J Math Biol. 1998, 37: 235252. 10.1007/s002850050127.
 11.
Roy P, Kozma R, Majumder D: From neurocomputation to immunocomputation: Model and algorithm for fluctuationinduced instability in biological systems. IEEE Trans. Evol. Computation. 2002, 9 (3): 292305.
 12.
Bizzari M, Cucina A, D’Anselmi F: Beyond the oncogene paradigm: understanding complexity in cancerogenesis. Acta Biotheor. 2008, 56: 173196. 10.1007/s1044100890478.
 13.
Kafarov V: Cybernetic methods in chemistry and chemical engineering. 1986, Moscow: Science Publishers
 14.
Singh S, Padhi R: Automatic path planning and control design for UAV using dynamic inversion. Proc. Amer. Control Conf. 2009, St. Louis: American Automatic Control Council, 24092414.
 15.
Khalil H: Nonlinear Control Systems. 1996, Prentice Hall: Princeton, 1996
 16.
Padhi R, Kothari M: An optimal dynamic inversionbased neuroadaptive approach for treatment of chronic myelogenous leukemia. Comp Methods and Prog in Biomed. 2007, 87 (3): 208224. 10.1016/j.cmpb.2007.05.011.
 17.
Padhi R, Balakrishnan S: Development and analysis of a feedback treatment strategy for parturient paresis of cows. IEEE Trans. Control Systems Tech. 2004, 12 (1): 5264. 10.1109/TCST.2003.821962.
 18.
de Pillis L, Gu W, Radunskaya A: Mixed immunotherapy and chemotherapy of tumor: Modeling, application and biological interpretations. J Theor Biol. 2006, 238: 841862. 10.1016/j.jtbi.2005.06.037.
 19.
de Pillis L, Radunskaya A, Wiseman C: A validated mathematical model of cellmediated immune response to tumor growth. Cancer Res. 2005, 65 (17): 79507958.
 20.
Calabresi P, Schein P: Medical Oncology: Basic Principles and Clinical Management of Cancer. 1993, New York: McGrawHill, 2
 21.
Yates A, Callard R: Cell death and the maintenance of immunological memory. Discrete & Contininous Dynamical Systems. 2002, 1: 4359.
 22.
Lanzavecchia A, Sallusto F: Dynamics of Tlymphocyte responses: intermediates, effectors, and memory cells. Science. 2000, 290: 9297. 10.1126/science.290.5489.92.
 23.
Stein A, Demuth T, Mobley D, et al: A mathematical model of gliobalstoma tumor spheroid invasion in a 3D in vitro experiment. Biophys J. 2007, 92: 356365. 10.1529/biophysj.106.093468.
 24.
Diefenbach A, Jensen E, Jamieson A, Raulet D: The Rael and H60 ligands of NKG2D receptor simulate tumor immunity. Nature. 2001, 413: 165171. 10.1038/35093109.
 25.
Dudley M, Wunderlich J, Robbins P, et al: Cancer regression and autoimmunity in patients after clonal repopulation with antitumor lymphocytes. Science. 2002, 298: 850854. 10.1126/science.1076514.
 26.
Macintyre E, Linch D: Lymphocytosis. Postgrad Medical J. 1988, 64: 4247. 10.1136/pgmj.64.747.42.
 27.
Jarosz M, Hak L, Wieckiewicz J, Mysliwska J: The NK cells in children with acute lymphoblastic leukemia and nonHodgkin lymphoma after intensive chemotherapy. Central Eur J Immunol. 2009, 34 (2): 9499.
 28.
Berrington J, Barge D, Spickett G: Lymphocyte subsets in term and significantly preterm UK infants. Clin Exp Immunol. 2005, 140 (2): 289292.445. 10.1111/j.13652249.2005.02767.x.
 29.
Jawahar S, Moody C, Chan M, Chatila T: Natural killer (NK) cell deficiency associated with an epitopedeficient Fc receptor type IIIA. Clin Exp Immunol. 1996, 103: 408413.
 30.
Roitt I, Brostoff J, Male D: Immunology. 1993, St. Louis: Mosby
 31.
Zotin A: Thermodynamic Bases of Biological Processes. 1990, New York: Walter de Gruyter Press
 32.
Ueno T, Ko S, Grubbs E, et al: Modulation of chemotherapy resistance in regional therapy: a novel therapeutic approach to advanced extremity melanoma using intraarterial temozolomide in combination with systemic O6benzylguanine. Mol Cancer Ther. 2000, 5: 732738.
 33.
Glick R, Lichtor T, Mogharbel A: Intracerebral vs. subcutaneous immunization with allogeneic fibroblasts genetically engineered to secrete interleukin2 in the treatment of glioma and melanoma. Neurosurgery. 1997, 41: 898907. 10.1097/0000612319971000000025.
 34.
Ganguly R, Puri I: Mathematical model for chemotherapeutic drug efficacy in arresting tumour growth based on the cancer stem cell hypothesis. Cell Prolif. 2000, 40: 338354.
 35.
Michor F, Hughes T, Iwasa Y: Dynamics of chronic myeloid leukaemia. Nature. 2005, 435: 12677031. 10.1038/nature03669.
 36.
Ewend M, Thompson R, Anderson R: Intracranial paracrine interleukin2 therapy stimulates prolonged antitumor immunity that extends outside CNS. J Immunotherapy. 2000, 23 (4): 438448. 10.1097/0000237120000700000007.
 37.
Kim TG, Kim CH, Park JS, Park SD, Kim CK, Chung DS, Hong YK: Immunological factors relating to the antitumor effect of temozolomide chemoimmunotherapy in a murine glioma model. Clin Vaccine Immunol. 2010, 17 (1): 143153. 10.1128/CVI.0029209.
 38.
SanchezPerez LA, Choi BD, Archer GE, Cui X, Flores C, et al: Myeloablative temozolomide enhances CD8^{+} Tcell responses to vaccine and is required for efficacy against brain tumors in mice. PlosOne. 2013, 8 (3): e5908210.1371/journal.pone.0059082.
 39.
Cheema TA, Wakimoto H, Fecci PE, Ning J, Kuroda T, Jayaretna DS, Martuza RL, Rabkin SD: Multifaceted oncolytic virus therapy for gliobalstoma in an immunocompetent cancer stem cell model. Proc Natl Acad Sci U S A. 2013, 110 (29): 1200612011. 10.1073/pnas.1307935110.
 40.
Bowles A, Perkins E: Longterm remission of malignant brain tumors. Neurosurgery. 1999, 44 (3): 636642. 10.1097/0000612319990300000110.
 41.
Lasky JL, Panosyan EH, Plant A, Davidson T, Yong WH, Prins RM, Liau LM, Moore TB: Autologous tumor lysatepulsed dendritic cell immunotherapy for pediatric patients with newly diagnosed or recurrent highgrade gliomas. Anticancer Res. 2013, 33 (5): 20472056.
 42.
Jones C, Schlosser M, van Zijl P: Amide proton transfer imaging of human brain tumors at 3 Tesla. Magnetic Reson Med. 2006, 56: 585592. 10.1002/mrm.20989.
 43.
Niranjan M, Liu X: State and parameter estimation of the heat shock response system using Kalman and particle filters. Bioinformatics. 2012, 28: 15011507. 10.1093/bioinformatics/bts161.
 44.
Baker S, Poskar C, Schreiber F, Junker B: An improved constraint filtering technique for inferring hidden states and parameters of a biological model. Bioinformatics. 2013, 29 (8): 10521059. 10.1093/bioinformatics/btt097.
 45.
Hidalgo M, Bruckheimer E, Rajeshkumar N: A pilot clinical study of treatment guided by personalized tumorgrafts in patients with advanced cancer. Mol Cancer Therapy. 2011, 10 (8): 13111316. 10.1158/15357163.MCT110233.
 46.
Disis M, Bernhard H, Jaffee E: Use of tumourresponsive Tcells as cancer treatment. Lancet. 2009, 373: 673683. 10.1016/S01406736(09)604049.
 47.
Monte U: Letter to the Editor. Cell Cycle. 2009, 8 (3): 505506. 10.4161/cc.8.3.7608.
 48.
Khar A: Mechanisms involved in natural killer cell mediated target cell death leading to spontaneous tumour regression. J Bioscience. 1997, 22 (1): 2331. 10.1007/BF02703615.
 49.
Lucas P, McNeil N, Hilgenfeld E, Gress R: Transforming growth factor β pathway serves as a primary tumor suppressor in CD8+ Tcell tumorigenesis. Cancer Res. 2004, 64: 65246529. 10.1158/00085472.CAN040896.
 50.
Baitsch L, Baumgaertner P, Devevre E, Speiser D: Exhaustion of tumorspecific CD8 + T cells in metastases from melanoma patients. J Clin Invest. 2011, 121 (6): 23502360. 10.1172/JCI46102.
Acknowledgements
The research was supported by core funds from National Brain Research Centre, Ministry of Science & Technology, Govt. of India. SK, VPSR and CS are respectively supported by funding received from educational grant of Dept. of Biotechnology, Dept. of Information Technology, and Ministry of Human Resource Development, Government of India.
Author information
Additional information
Competing interests
The authors declare that they do not have any competing interests.
Authors’ contributions
The research was conceived, designed and planned by PR and RP. The modelling and numerical experiments were performed and the biological studies investigated by SK, VPSR, CK, RP and PR. All authors examined and evaluated the data. SK, VPSR and PR wrote the manuscript, with participation from CS and RP. All authors read and approved the final manuscript.
Suhela Kapoor, VP Subramanyam Rallabandi contributed equally to this work.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Received
Accepted
Published
DOI
Keywords
 Tumour cell extinction
 Cancer therapy optimization
 Control system
 Chemotherapy
 Immunotherapy
 Glioma
 Astrocytoma
 Meningioma
 Oligodendroglioma
 Glioblastoma