On systems and control approaches to therapeutic gain

Simulation

Systems of first order nonlinear ordinary differential equations (ODEs) [10] are commonly used to model and simulate biochemical networks. The models have boundary condition nodes that decouple them from their environment, and state variable nodes that are responsive to environmental changes. The models are typically built and validated using 4 types of data:

1. Configuration data – These data specify how various system components are interconnected. Wiring diagrams are an excellent example of this form of data.

2. Isolated component data – These data describe how system components behave in complete isolation, e.g. enzyme mechanisms and estimates of binding constants.

3. Operating point data – These data describe the state of the interconnected system in steady state, e.g. metabolite concentrations in whole cell extracts of unperturbed cells.

4. Dynamic response data – These data specify how the interconnected system reacts to perturbations, e.g. metabolite concentration time course responses to drugs.

In general, configuration data, isolated component data and operating point data are used to build models, and dynamic response data are then used to validate them, e.g. see [3]. Once built, the models can be represented using Systems Biology Markup Language (SBML) [11] and analyzed using various software packages [12], such as R [13] and Matlab [14].

Control

The inverse of simulation is control. The goal in control is to identify system inputs (boundary condition time courses in this paper) that will drive the system toward an objective, such as a target steady-state (a regulation problem) or a desired time course trajectory (a tracking problem). Various types of control methods can potentially be applied to problems of interest in cancer research, e.g. classical techniques such as proportional, integral and derivative (PID) control, or model-based methods which involve state estimation as a component of the feedback system. In model-based methods, the first step is to develop a validated process model, and the second step is to use the model and the overall problem objective to design a control strategy. In the future, as cancer treatment protocols are eventually developed based on models and control system design methodologies (Fig. 1), and as model parameters and initial states become individualized through patient specific laboratory measurements, treatment dose time course schedules will become individualized as a direct consequence of using model-based control system designs. Further, since successful control systems in nature and in industry are generally implemented using closed-loop feedback strategies, where current state estimates derived from available observations are used to determine the most suitable input values in the next time interval, as science advances and as measurement costs decrease, it is extremely likely that state feedback based cancer treatments will eventually become the standard of care. This suggests that systems cancer biology should perhaps focus first where it is most likely to perform the best, namely, on problems where control system techniques can most readily be applied. In this regard, leukemias are attractive to study because leukemic cells can be sampled to measure relevant molecular activities [15–22] in real time (i.e. daily), as needed to implement feedback control strategies. Leukemias are also attractive to study because they seem to be simpler than solid tumors, perhaps because white blood cells naturally leave and enter tissues and thus require fewer steps to become malignant, a notion supported by shorter latency times for radiation induced leukemias [23–25] versus solid tumors [26]. Since increased process understanding favors process controllability, leukemia simplicity relative to solid tumors should favor its curability.

In the past, control system design methods have been developed primarily for linear [27] rather than nonlinear systems [28], and most often for systems with rich rather than sparse time course data. New design methods and/or ad hoc solutions will therefore be needed to successfully apply control system approaches to cancer therapy.

Previous attempts to apply control theory [29, 30] to cancer therapy [31, 32] lacked sufficient data, and the process models that were used were not molecularly based. As our understanding of cancer relevant biological processes increases, our knowledge of the dynamic behavior of these processes will improve, and thus so too will the likelihood of successfully applying control system design approaches to cancer therapy. The purpose of this paper is to formulate (but not solve) two systems and control oriented, biochemical system model-based therapeutic gain strategies as abstractions of two conceptually distinct, specific cancer treatment strategies.

Direct approach to therapeutic gain

Iododeoxyuridine (IUdR) sensitization of cells to radiation induced cell killing [39, 40] is proportional to its incorporation into DNA [41, 42], and IU-DNA has a longer half-life in mismatch repair deficient (MMR^-) cells than in MMR⁺ cells [43, 44]. Thus, IU-DNA levels at the time of irradiation correlate with the probability of clonogenic cell death and are modulated by the cause of the cancer. Suppose we have a model of IUdR metabolism that includes IU-DNA and MMR. If the MMR parameters of this model are set to MMR deficient values, the model represents MMR^- malignant cells, otherwise, with wild type MMR parameters, the model represents MMR⁺ normal cells. Such a model, or pair of models, could be used to develop novel multi-drug approaches to the treatment of MMR^- cancers. Specifically, the models could be used to predict how IUdR metabolism should be manipulated by drugs over time to maximize IU-DNA differences between MMR^- and MMR⁺ cells at some time point. The timing of the maximum difference would then determine the best timing of acute irradiation to exploit the IU-DNA differences for a therapeutic gain. This treatment scenario is depicted in Figure 2.

The problem statement outlined above is direct to the extent that its goal is selective killing of malignant versus normal cells. It requires a model that relates anti-cancer agents to an endpoint that is modulated by the cancer causing event and that correlates with cell death (Fig. 3A). These requirements define a class of cancer therapy problem statements with a common abstraction. The abstraction involves a pair of tissue specific dynamical system models of the same biological process, both manipulated synchronously through a set of common input functions. Subject to normal tissue toxicity constraints, the overall objective is to maximize the expectation of differential cell killing, defined as the probability of killing a malignant cell relative to the probability of killing a normal cell (i.e. the odds), or the log thereof. This defines a class of optimal control problems and a systems and control approach to therapeutic gain referred to hereafter as the direct approach.

To begin to formalize the direct approach mathematically, consider the pair of system models $\dot{x}$ ⁺ (t) = f(x ⁺(t), u(t), p ⁺, $x_0^{+}$) and $\dot{x}$ ^- (t) = f(x ^-(t), u(t), p ^-, $x_0^{-}$) for normal and malignant cells, respectively, where x(t) is a vector of state variables (e.g. metabolite and protein concentrations), u(t) is a vector of input time courses applied synchronously to both systems (e.g. extracellular drug concentrations), p is the set of model parameters that differ between the two cell types (parameters with identical values in both cell types are embedded in f), and the initial states of the normal and malignant cells are $x_0^{+}$ and $x_0^{-}$, respectively. Here, normal cells may be dose-limiting marrow or gut cells, or they may be normal tissue counterparts of malignant cells. Since the two models are of the same biochemical system, they represent the same chemical species and thus have identical state space dimensions. Some of these dimensions will be more controllable than others, some will correlate with the probability of cell death more than others, and some, particularly those associated with major parameter differences, will separate values set initially to identical values in the two models more than others. For direct approach success, cancers of focus should be selected such that controllability, observability and separability are aligned within the same dimension(s) as much as possible.

Indirect approach to therapeutic gain

Pro-B cell childhood acute lymphoblastic leukemia (ALL) [45] has event-free five year survival probabilities of 85–90% for the TEL-AML subtype and 20–40% for the BCR-ABL subtype [46]. Thus, TEL-AML1 ALL is essentially curable, BCR-ABL ALL is essentially non-curable, and both have pro-B cell origins. Suppose a biochemical system model existed that connected drug targets to a critical (see below) determinant of BCR-ABL treatment failure. Given such a model, a treatment failure risk state transfer strategy can be envisioned where the goal is to identify drug concentration time course schedules that transfer non-curable BCR-ABL initial states toward a target state region defined by the collection of curable TEL-AML1 states. In this scenario, if poor prognosis BCR-ABL patients could be brought to a state in the neighborhood of curable TEL-AML1 states, they might then be treated effectively from such a region by a standard TEL-AML1 therapeutic regimen.

Because effective combination TEL-AML1 therapies include methotrexate (MTX), it was hypothesized that differences in folate system states may exist between TEL-AML1 and BCR-ABL patients. To explore this hypothesis, the diagnostic bone marrow data of Yeoh et al [36] was used in conjunction with the folate model of Morrison and Allegra [3] to develop model-based steady state flux predictions of de novo purine synthesis (DNPS) versus de novo thymidylate synthesis (DNTS) as shown in Figure 4. In this plot, if the center of the TEL-AML1 patient states serves as the target of BCR-ABL patient state space trajectories, at first sight, it appears that BCR-ABL patients might benefit by replacements of MTX with ALIMTA, since ALIMTA inhibits glycinamide ribonucleotide formyltransferase and thus DNPS more so than DNTS [47]. This illustration of the concept of the indirect approach (Figure 3B) must, however, be carefully scrutinized before it is used to suggest new treatments. For example, if the one cured BCR-ABL patient (B in Fig. 4) serves as the guide, the other BCR-ABL patients should then be driven toward higher DNPS and higher DNTS (upper right in Fig. 4), rather than lower DNPS and higher DNTS (lower right in Fig. 4). Further, the TEL-AML1 patients with hematological relapses (T in Fig. 4) do not fall within the BCR-ABL cluster, and several TEL-AML1 patients who do fall within the BCR-ABL cluster were cured. Thus, closing differences in this plot may have no impact on treatment outcome. This suggests that the critical determinants of differences in treatment outcome lie downstream of the folate system, or are associated with other drugs in the TEL-AML1 regimen.

BCR-ABL is a prognostic determinant of treatment outcome because it defines the poor prognosis group. BCR-ABL is also likely to be the cause of BCR-ABL leukemias. However, this does not imply that successfully closing BCR-ABL activity differences cures patients. Above and beyond being prognostic, a critical determinant of treatment outcome can be defined as one which leads to improved treatment outcomes after difference closure. For example, because Gleevec does not cure BCR-ABL ALL, BCR-ABL is not a critical determinant of BCR-ABL ALL. However, since ABL expression (which is equivalent to BCR-ABL expression because probe sets are taken from the 3' end) is minimal for the cured BCR-ABL patient in 3 out of 4 probesets with P < .0.06 differences between leukemia subtypes, Fig. 5 does suggest (very mildly) that Gleevec combined with standard therapy could be beneficial, i.e. that BCR-ABL may be conditionally or partially critical. Thus, efforts to model successful drugs such as MTX in TEL-AML1 ALL might be productively complemented by parallel efforts to model candidate critical determinants of outcome.

The indirect approach is a generalization/abstraction of the BCR-ABL example given above. It can be characterized as an initial state/target state control problem, formulated within the framework of a dynamical system model. This is reminiscent of the problem of ballistic controllability that forms the foundation of the notions of controllability and reachability in modern systems theory. The ballistic control problem is one in which the objective is to find an input signal (i.e. control) that transfers a given initial state of the system to a specified final state in a given period of time. The basic idea of ballistic controllability is best understood for the simple system defined by the linear system model: $\frac{dx\left(t\right)}{dt}$ = B(t)u(t) where x(t) is the state of the system (x ∈ Rⁿ), u(t) is the input (u ∈ R^m, with m < n) and B(t) is a time varying system input matrix. The objective is to transfer the given initial state x(t₀) = x₀ to the specified final state x(t_f) = x_f over the time interval [t₀, t_f], i.e. the problem is to determine the m time courses in u such that $x\left(t_f\right)={\displaystyle \underset{t_0}{\overset{t_f}{\int }}B\left(\tau \right)u\left(\tau \right)d\tau +x\left(t_0\right)}$. For this problem, we introduce the controllability gramian $W\left(t_0,t_f\right)={\displaystyle \underset{t_0}{\overset{t_f}{\int }}B\left(\tau \right)B^{\prime }\left(\tau \right)d\tau }$. If W(t ₀, t _f ) is full rank, then there exists an input signal vector u(t) that transfers x(t₀) = x₀ to x(t_f) = x_f . In fact, u(t) = B'(t)η where η = W(t₀, t_f)^-1(x_f-x₀), i.e. the desired state transfer input signal is determined uniquely by the parameters of the problem. This simple problem can be easily extended to the class of linear systems: $\dot{x}$(t) = A(t)x(t) + B(t)u(t)[27] and to the general class of nonlinear systems that are linear in the control effort u(t): $\dot{x}$(t) = f (x(t)) + g(x(t))u(t) [28]. Important questions remain. For example, should the system state be clamped to lie within the target state region while applying standard therapy? Or should it be released to avoid interference between clamp control efforts and standard therapy control efforts? Or should a compromise between these two extremes be sought?