Outcome Prediction in Mathematical Models of Immune Response to Infection

Clinicians need to predict patient outcomes with high accuracy as early as possible after disease inception. In this manuscript, we show that patient-to-patient variability sets a fundamental limit on outcome prediction accuracy for a general class of mathematical models for the immune response to infection. However, accuracy can be increased at the expense of delayed prognosis. We investigate several systems of ordinary differential equations (ODEs) that model the host immune response to a pathogen load. Advantages of systems of ODEs for investigating the immune response to infection include the ability to collect data on large numbers of ‘virtual patients’, each with a given set of model parameters, and obtain many time points during the course of the infection. We implement patient-to-patient variability v in the ODE models by randomly selecting the model parameters from distributions with coefficients of variation v that are centered on physiological values. We use logistic regression with one-versus-all classification to predict the discrete steady-state outcomes of the system. We find that the prediction algorithm achieves near 100% accuracy for v = 0, and the accuracy decreases with increasing v for all ODE models studied. The fact that multiple steady-state outcomes can be obtained for a given initial condition, i.e. the basins of attraction overlap in the space of initial conditions, limits the prediction accuracy for v > 0. Increasing the elapsed time of the variables used to train and test the classifier, increases the prediction accuracy, while adding explicit external noise to the ODE models decreases the prediction accuracy. Our results quantify the competition between early prognosis and high prediction accuracy that is frequently encountered by clinicians.


Introduction
The immune response to infection is a complex process that involves a wide range of length scales from proteins to cells [1][2][3][4], tissues [5], and organ systems [6].Despite enormous progress over the past 30 years in developing mathematical models for the immune response to infectious disease such as tuberculosis [7], HIV [8][9][10][11], and influenza [12,13], these models still have not been able to dramatically improve patient diagnosis, prognosis, and treatment [14,15].Instead, vaccine and drug development often relies on costly trial-and-error methods [16].However, advances in gene sequencing capabilities [17], increasing speeds of computer processors, and the ability to store enormous amounts of medical data promise dramatic improvements in mathematical approaches to predicting and controlling the response to infectious disease [18][19][20].
One promising mathematical approach is to use machine learning methods on large data sets to classify patients as healthy or sick, perform early warning analyses for early detection of infection, or identify the minimal set of genes responsible for a particular immune response.[21,22] However, many questions are left unanswered in such studies.For example, how much and what kinds of data are required to have confidence in the machine learning predictions and what are the underlying biophysical mechanisms for the relationships between variables that are identified by these techniques?Further, it is difficult to determine differences in the immune response that arise from patient-to-patient variations compared to slight differences in the initial conditions of each patient.
In this manuscript, we focus on sets of ordinary differential equations (ODEs) as mathematical models for the immune response to infection.The advantages of ODEs are manifold: 1) Each 'virtual patient' can be considered as a set of parameters in the set of ODEs; 2) There is essentially no limit on the amount of data that can be collected on each virtual patient; 3) The accuracy of machine learning predictions can be explicitly tested as a function of the number of time points and initial conditions for each patient and the number of patients included in the training and testing sets; and 4) analysis of the fixed points (or steady-state outcomes) and basins of attraction of the ODEs can give biophysical insight into the immune response to infection.
We will investigate several classes of ODE models for the immune response to infection.First, we will describe a four-dimensional model for the acute inflammatory response to a pathogen load that was studied in detail in Ref. [23].We will then consider reduced versions of this model with fewer variables and parameters obtained by slaving one or more of the original four variables, as well as changes to form of the ODEs that alter the fixed point structure and flows between them.For each model, a virtual patient is defined by one set of parameters.Given an initial condition (values of the variables at time t = 0), the patient will evolve deterministically to one of several possible discrete steady-state (t → ∞) outcomes, or fixed points.Thus, for each patient, we can determine the basins of attraction that map initial conditions for all of the variables to steady-state outcomes by numerically integrating the sets of ODEs.
We seek to determine the limits of the prediction accuracy of discrete steady-state outcomes of ODEs as a function of patient variability (i.e.random fluctuations in parameter values) using machine learning techniques.In the limit of zero patient variability, our simple classification algorithm (logistic regression) can achieve nearly perfect prediction accuracy even when the classification occurs on variables at short times.However, as the patient variability increases, the basins of attraction for different patients yield different outcomes for a given initial condition as shown in Fig. 1 for 2% patient variability in model (1) for the immune response to infection.(See Materials and Methods.)The fact that each set of initial conditions does not possess a unique outcome places a fundamental limit on the predictability of patient outcomes.Thus, we find that the machine learning prediction accuracy decreases with increasing patient variability.In contrast, for a given patient variability, the prediction accuracy increases with the time used for classification as the systems converge to their steady-state outcomes (Fig. 1).We also show that at short times our classification algorithm saturates the theoretical limit for the prediction accuracy in the presence of patient variation, and that the addition of external noise only worsens the outcome prediction accuracy.
The manuscript is organized as follows.In the Materials and Methods section, we introduce several ODE systems that have been used to model the host immune response to infection [23], including their parameter sets and discrete steady-state outcomes, and describe how we implement patient variability in the ODE models.In the Results section, we emphasize our three main results that hold for all of the ODE models we studied: 1) patient variability leads to overlap of the basins of attraction for the steady-state outcomes, which limits the outcome prediction accuracy, 2) the prediction accuracy increases with the time used for classification because the basins of attraction separate with increasing time, and 3) the addition of external measurement noise further reduces the prediction accuracy.In the Discussion section, we point out the clinical implications of our work and describe important future studies of the prediction accuracy for ODE models with continuous outcomes.In the Supporting Information, we discuss a generalized one-dimensional ODE against which we compare our results and the numerical implementation of the algorithms used to classify the steady-state outcomes of the ODEs.

Materials and Methods
Our studies focus on several ODE models with varying complexity for the immune response to pathogen load that were first introduced in Ref. [23].These ODE models can have up to four coupled variables that represent the concentration of pathogen P , activated neutrophils N , inflammation (or damage) D, and immuno-suppressor (cortisol) C. The models include interactions between these four quanties.For example, the presence of pathogen P > 0 causes an immune response, where neutrophils are activated and N increases.Neutrophils kill pathogen, which decreases P , but also cause inflammation (damage), which increases D. The cortisol level C increases when there is a high neutrophil level, which then reduces the neutrophil level.
Model (1) (Eqs.1-4) includes all four variables P , N , D, and C. The right-hand side of dP/dt is a sum of three terms.The first term enables logistic growth of the pathogen.In the absence of any other terms, any positive initial P 0 will cause P to grow logistically to the steady-state value P ∞ .The second term mimics a local, non-specific response to an infection.For small values of P , the decrease is proportional to P .For larger values of P , the decrease caused by the second term is constant.The third term models the decrease of P due to interactions with activated immune cells (neutrophils) N .Activated neutrophils N can directly decrease P .The anti-inflammatory response, which is captured by the cortisol level C, mitigates this effect leading to a decrease of P proportional to N * P/(1 Two terms determine the rate of change in neutrophils, dN/dt.The first term accounts for the fact that neutrophils can be activated if a resting neutrophil cell encounters a pathogen P or an already activated neutrophil N .Furthermore, tissue damage D also triggers the activation of neutrophils.The second term describes the death of neutrophils N , with the decrease in N proportional to the amount of neutrophils present. The rate of change in damage dD/dt is also controlled by two terms.The first term mimics positive feedback between D and N .Activated phagocytes cause collateral damage in the tissue.Again, the effectiveness of N is mitigated by the anti-inflammatory response 1/(1 + (C/C ∞ ) 2 ).The saturation function f s models the fact that the effect of N on D saturates for large N .The second term, -µ d D, represents repair of the tissue.
The anti-inflammatory response C increases with the source term s c .In addition, there is a natural death rate µ c , which leads to a positive steady-state value of C in the absence of any immune activation N or damage D. However, even small amounts of damage and neutrophils will up-regulate C. In the case of small N + k cnd D, the production of C is proportional to N + k cnd D, while for large values of N + k cnd D, changes in C are proportional to k cn .Again, the effectiveness of N is mitigated by Model (1) has 21 parameters: and µ c .Depending on the values of these parameters, model (1) possesses different numbers of fixed points with varying stabilities.However, we will focus on a specific parameter regime (given in Table 1) with three stable fixed points, which correspond to the physiological steady-state outcomes: health, septic death, and aseptic death.

Model (1)
where  2)-( 5) given below are simplified versions of model (1).A summary of the dimension, number of parameters, and number of stable fixed points for each of the ODE models is shown in Table 2. To obtain model ( 2) from (1), C is set to a constant C = 0.23 and the remaining terms define a three-variable model with P , N , and D. For model (3), we set C = C and D = 0, which gives a two-variable model for P and N .For model (4), we set C = C and P = 0 to obtain a two-variable model for N and D. In this model, the value of the initial rise in N can be thought of as the response to trauma.For model (5), we set C = C = 0.1, D = 0, and N = 0, which gives a one-dimensional model for P .This model only treats the innate immune response with no activated neutrophils.
Model ( 3) where Model ( 4) where Model ( 5) In Fig. 2, we show the time evolution of the four variables P , N , D, and C for model (1) for twenty different sets of random initial conditions to illustrate its three stable fixed points (health, septic death, and aseptic death) using the parameter values in Table 1.For trajectories that approach the septic death fixed point, the pathogen and neutrophil levels grow rapidly.The high neutrophil level causes cortisol to increase as well.Despite the high level, the neutrophils cannot reduce the pathogen load and the cortisol level is not large enough to reduce the neutrophil level.As a result, the high neutrophil level causes significant damage at long times, which is termed septic death due to the associated high pathogen level.Thus, the septic death steady-state outcome is characterized by P > 0, N > 0, D > 0, and C > C ∞ .
In the healthy state, the pathogen level can be reduced to zero by the neutrophils, and the neutrophil level can be reduced to zero by cortisol.Once the neutrophil level is zero, the cortisol level returns to its background level and damage decreases to zero.Thus, the healthy state is characterized by P = 0, N = 0, D = 0, and C = C ∞ .For the parameter values in Table 1, model (1) possesses three fixed points: health (green lines), septic death (purple lines), and aseptic death (orange lines).The initial conditions are sampled randomly within the cube: 0 ≤ P 0 ≤ 0.42, 0 ≤ N 0 ≤ 0.255, D 0 = 0, and 0 ≤ C 0 ≤ 0.35.The three fixed points can be differentiated by the steady-state values of P and D: health (P = 0, D = 0), aseptic death (P = 0, D > 0), and septic death (P > 0, D > 0). 10, 7, and 3 of the initial conditions evolve to the health, septic death, and aseptic death fixed points, respectively.1) with 20 sets of randomly selected parameters with the same initial conditions.10% patient variation allows the system to reach the health and septic death fixed points with the initial condition P 0 = 0.45, N 0 = 0.45, D 0 = 0, and C 0 = 0.35, whereas only the aseptic death fixed point is obtained for this initial condition with no patient variation.5, 4, and 11 of the trajectories evolve to the health, septic death, and aseptic death fixed points, respectively.
During the approach to the aseptic death fixed point, the neutrophil level is strong enough to reduce the pathogen level to zero, but the cortisol level is insufficient to reduce the neutrophil level to zero, which leads to increasing damage.Thus, the aseptic death fixed point is characterized by P = 0, N > 0, D > 0, and C > C ∞ .
The outcomes of the immune response to infection can vary from patient to patient, even with the same initial conditions (e.g. the pathogen load).To introduce patient variability into the ODE models, we select the parameters ({q}) in models (1)-( 5) randomly from independent Gaussian distributions with mean values µ q in Table 1 and variance v relative to the mean.Negative values of the parameters can cause the ODE models to become non-integrable, and thus the parameter distributions are cut off so that the parameter values are non-negative.We solve the ODE models for 10 4 sets of parameters for each of the 10 4 random initial conditions at each v.The limits for the sampling of the initial conditions for each model are given in Table 3.We then perform a classification analysis on these trajectories to predict the steady-state outcomes.The prediction accuracy A is defined as the number of correct classifications of the steady-state outcomes divided by the total number of classifications.

Results
For a deterministic system of ODEs, the basin of attraction for a given fixed point is defined as the collection of initial conditions that evolve to that particular fixed point.For a given set of parameters, each of the ODE models ( 1)-( 5) possesses well-defined (non-overlapping) basins of attraction for each fixed point.However, different outcomes can be achieved even for a single initial condition if the parameters of the ODE model are varied.(See Fig. 3.) For example, the ratio of the parameters s c and µ c determines the background level of cortisol in model (1).Background cortisol levels are known to vary from patient to patient and can vary from one organ system to another in a given patient.To mimic these variations, we select sets of parameters randomly with mean values in Table 1  The prediction accuracy A using a logistic regression classifier at time t c = 0 (symbols) and the average best guess over 10 3 initial conditions (dashed curves) versus patient variation v for (a) models ( 1) and ( 2), (b) models ( 3) and (4), and (c) models ( 5) and (6).
We seek to predict the patient steady-state outcomes in models (1)-( 5) in the presence of patient variability v.For the prediction method, we employ logistic regression with one-versus-all classification [24].We compare the prediction accuracy at patient variability v to the average best guess of the steady-state outcome.For the best guess method, we determine the steady-state outcome for each of 10 2 sets of parameters for a given initial condition.We define the best guess as the steady-state outcome with the highest number of occurrences and record the frequency f i of the best guess for initial condition i.We then average the frequency f i over 10 3 initial conditions for each v to obtain an estimate for the prediction accuracy in systems with basin overlap.
For the prediction method, we solve a given system of ODEs for N i = 10 4 random initial conditions, each with randomly selected parameter sets with variance v.We choose N t = 800 of the N i trajectories randomly to train the classifier and predict the outcome of the remaining 9200 trajectories.The classifier maps the state of the system at a given time t c to a particular steady-state outcome.The prediction accuracy is then averaged over 10 training and prediction runs, each with N t = 800 randomly selected training trajectories.1) and ( 2), (b) models ( 3) and ( 4), and (c) models ( 5) and ( 6) for patient variation v = 0.05.
In Fig. 4, we compare the accuracy A of the logistic regression prediction method (with classification at time t c = 0) to the average best-guess frequency as a function of the patient variability v for models (1)-( 5).For all model ODEs, the prediction accuracy for the logistic regression prediction method is near 100% at v = 0, decreases for increasing patient variability, and reaches a plateau near 1/n f in the large v limit, where n f is the number of stable fixed points in the model (except for model ( 3)).For model (3) with two steady-state outcomes, the prediction accuracy is non-monotonic and increases for v > 0.3 because this ODE model begins to sample parameter regimes where one steady-state outcome is much more probable than the other.In addition, for all models the average best-guess frequency provides an upper bound for the accuracy of the prediction algorithm.Hence, the overlap of the basins of attraction imposes a limit on the prediction accuracy.To test the generality of these results, we studied another one-dimensional ODE (model ( 6)) with varied fixed point structure compared to that for model (5).(See Supporting Information.)The results for model (6) are very similar to those for models (1)- (5).
In Fig. 4, we showed results for the logistic regression prediction method with classification at t c = 0.In Fig. 5, we show the prediction accuracy for models (1)- (6) with patient variability v = 0.05 as a function of the classification time t c .For all models, the prediction accuracy grows with increasing t c , reaching nearly 100% beyond a characteristic time t * that depends on the model.The prediction accuracy improves at later classification times because the system trajectories have moved closer to the fixed points and hence the basins of attraction are more easily separated as shown in Fig. 1 for model (1).We also investigated the variation of the prediction accuracy in the presence of measurement noise.We took the trajectories generated for Fig. 4 and added Gaussian random noise to the model variables with variance s at each time point.We then performed training and testing on the noisy data with classification at time t c = 0.In Fig. 6, we show for model ( 6) that the prediction accuracy decreases with increasing s.We find similar results for models (1)-( 5).These results emphasize that even if the measurement noise could be reduced to zero, the patient variation imposes an intrinsic limitation to outcome prediction.

Discussion
In clinical settings it is of great importance to determine patient outcomes as quickly as possible with maximum accuracy.In this manuscript, we studied the effects of patient variability on the ability to predict steady-state outcomes in systems of ODEs that model the immune response to infection.For deterministic systems of ODEs with a given fixed set of parameters, each initial condition can be mapped to a given steady-state outcome (or fixed point) and the collection of initial conditions that map to a given steady-state outcome is defined as the basin of attraction of that outcome.Each virtual patient can be defined by a given set of parameters in the model ODE and patient variability can be introduced by varying the model parameters.
We showed that the introduction of patient variation leads to overlaps of the basins of attraction for the steady-state outcomes.In particular, a given initial condition can map to multiple steady-state outcomes for different virtual patients (i.e.v > 0), which is similar to the case of patients showing different responses to infection in clinical settings.We find that the prediction accuracy of the outcomes decreases strongly with increasing patient variability.Our results emphasize that even when the complete state of the system is known (i.e.all patient variables are measured precisely as a function of time), we have limited knowledge of the patient outcome when there is patient-to-patient variability that gives rise to basin overlap.
Our results also show that for all of the model ODEs studied the prediction accuracy increases as the time t c used for classification increases.As t c increases, the systems move closer to their steady-state outcomes and the basins of attraction separate, which increases the prediction accuracy.Again, this result is consistent with clinical experience.If a clinician waits to see if the condition of the patient improves or worsens, the prognosis will become more accurate.In our work, we explicitly show that patient-to-patient fluctuations cause a competition between early and accurate outcome prediction.
In this work, we focused on discrete steady-state outcomes (i.e.health or death of the patient) of the immune response to infection.However, in many biomedical scenarios, the outcomes involve continuous variables rather than discrete states.In future work, we will apply similar techniques to understand the effects of patient variability on the predictions of continuous model variables, for example, the immune response and vaccination efficacy for influenza [25].
we label the outcomes, y (i) = 0 and y (i) = 1 for the ith set of variables.We fit a logistic (or sigmoidal) function to the input data such that the function gives the probability that y (i) = 1 given the input data x (i) .The parameters θ j , where j = 1, 2, . . ., N , are determined by minimizing the cost function J(θ 1 , θ 2 , . . ., θ N ) = − 1 m m i=1 y (i) log(P θ (x (i) )) + (1 − y (i) )log(1 − P θ (x (i) )) , where m is the number of training samples.Since models (1) and ( 2) possess three steady-state outcomes (aseptic death, septic death, and health), we must go beyond the binary classification scheme described above.To classify ODE models with three steady-state outcomes, we implement the one-versus-all classification scheme.To do this, we consider three labeled outcomes, y h , y ad , and y sd , for a given set of variables x (i) .y (i) h = 0 if the patient outcome is not health (i.e.aseptic or septic death) and y (See Table 4.) We use these labeled outcomes and Eq.18 to determine P h (x), P ad (x), and P sd (x).Given an unlabeled set of variables (x = {x 1 , x 2 , . . ., x N }), we calculate P h , P ad , and P sd and select the outcome with highest probability to be the predicted outcome.

Figure 1 .
Figure 1.Time evolution of patient outcomes for a range of neutrophil (N ) and cortisol (C) initial conditions for Model A. (left) Patient outcomes in the long-time limit given N and C initial conditions for model (1) are shaded green, orange, and purple for the steady-state outcomes of health, aseptic, and septic death, respectively.The initial values of the pathogen load and damage are P 0 = 0.35 and D 0 = 0, and patient variability is set to v = 2%.The right six panels indicate how the systems in the leftmost panel separate in the N and C plane as time increases, t = 10, 20, 50, 100, 250, and 500, from left to right.

Figure 3 .
Figure 3. P , N , D, and C versus time for model(1) with 20 sets of randomly selected parameters with the same initial conditions.10% patient variation allows the system to reach the health and septic death fixed points with the initial condition P 0 = 0.45, N 0 = 0.45, D 0 = 0, and C 0 = 0.35, whereas only the aseptic death fixed point is obtained for this initial condition with no patient variation.5, 4, and 11 of the trajectories evolve to the health, septic death, and aseptic death fixed points, respectively.

Figure 4 .
Figure 4. Prediction accuracy of the steady-state outcome as a function of patient variation.The prediction accuracy A using a logistic regression classifier at time t c = 0 (symbols) and the average best guess over 10 3 initial conditions (dashed curves) versus patient variation v for (a) models (1) and (2), (b) models (3) and (4), and (c) models (5) and(6).

Figure 5 .
Figure 5. Prediction accuracy of the steady-state outcome as a function of the classification time.The prediction accuracy A using a logistic regression classifier at time t c (symbols) for (a) models (1) and (2), (b) models (3) and (4), and (c) models (5) and (6) for patient variation v = 0.05.

Figure 6 .
Figure 6.Prediction accuracy as a function of patient variation for different noise strengths.The prediction accuracy A using a logistic regression classifier at time t c = 0 for model(6) in the presence of measurement noise with strength s = 0 (circles), 0.05 (squares), 0.10 (diamonds), 0.20 (triangles), and 0.50 (triangles).
(i) h = 1 if the patient outcome is health.Similar definitions apply for y (i) ad and y (i) sd .

Table 2 .
Summary of ODE modelsModel Dimension Parameters q Stable Fixed Points k pg