 Research
 Open Access
 Published:
Uncovering the effects of heterogeneity and parameter sensitivity on withinhost dynamics of disease: malaria as a case study
BMC Bioinformatics volume 22, Article number: 384 (2021)
Abstract
Background
The fidelity and reliability of disease model predictions depend on accurate and precise descriptions of processes and determination of parameters. Various models exist to describe withinhost dynamics during malaria infection but there is a shortage of clinical data that can be used to quantitatively validate them and establish confidence in their predictions. In addition, model parameters often contain a degree of uncertainty and show variations between individuals, potentially undermining the reliability of model predictions. In this study models were reproduced and analysed by means of robustness, uncertainty, local sensitivity and local sensitivity robustness analysis to establish confidence in their predictions.
Results
Components of the immune system are responsible for the most uncertainty in model outputs, while disease associated variables showed the greatest sensitivity for these components. All models showed a comparable degree of robustness but displayed different ranges in their predictions. In these different ranges, sensitivities were wellpreserved in three of the four models.
Conclusion
Analyses of the effects of parameter variations in models can provide a comparative tool for the evaluation of model predictions. In addition, it can assist in uncovering model weak points and, in the case of disease models, be used to identify possible points for therapeutic intervention.
Background
Malaria is a wellknown parasitic disease caused by Plasmodium parasites and affects populations in tropical and subtropical areas with a significant impact in the subSaharan African region. Although an approximate of 228 million cases of infection and a mortality rate of 405,000 has been reported for 2018, the World Health Organization states that the overall incidence of malaria has decreased from 2010 [1]. However, this rate of decrease plateaued with the development of resistance to current treatments, necessitating new investigative approaches into disease treatment and eradication [1, 2]. A wellstudied part of the parasite lifecycle is the blood stage within the human host [3]. In this erythrocytic stage, merozoites infect red blood cells, where they mature and proliferate until the red blood cells burst and release more merozoites [4]. The erythrocytic stage is associated with clinical symptoms of malaria since parasites and infected red blood cells (iRBCs) activate the immune system’s response, which can lead to symptoms such as fever, malaise and exacerbate anaemia [5]. The immune system, however, plays an important role in disease progression and response to treatment and vaccination. Taking its response into consideration is therefore vital for a quantitative understanding of infection dynamics and intervention effects [5].
Upon malaria infection, the alwayspresent innate immune system rapidly attempts elimination of the invading pathogen [6, 7]. The innate immune system is nonspecific to an infection and fights microbes by secreting various proteins and cytokines and assisting the adaptive immune system [8]. The adaptive immune response takes longer to develop and can be subdivided into two categories, namely humoral and cellmediated immune responses [7,8,9]. Here the humoral response would target free pathogens in the blood using antibodies produced from Blymphocytes, whereas the cellmediated response targets infected cells with Tcells differentiated from Tlymphocytes. Even with these multiple immune responses, malaria infection can persist in the body leading to severe complications and death [10].
Experimental methods such as bioluminescent, behavioural studies and a variety of assays can be used to investigate the malaria parasites, their interactions with the human and mosquito hosts and the immune system’s response to infection [6, 11, 12]. Clinical studies are also used to investigate parasite interactions and treatment effects, using methods such as microscopy, rapid diagnostic tests and molecular assays [13]. Data from clinical and experimental observations can be encoded in mathematical models. These models can then be used to analyse the behaviour emerging as a function of interacting processes in the system, and to make quantitative predictions of how a system is influenced by various alterations. A variety of withinhost mathematical models describing the disease dynamics associated with malaria infection exist, with differential equation based models focusing on the time evolution of different cell populations within the host [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]. Simple models of the erythrocytic infection stage usually describe populations of healthy red blood cells (RBCs), infected red blood cells (iRBCs) and free roaming merozoites with some of these models extended to include immune system components [21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37].
When comparing these models it becomes evident that they often differ in their formalism, structure and in the interpretation of model components and simulation results depending on the purpose of the original study. In each model, the parametrisation of the biological processes approximate the dynamics that might only be realistic close to the reference state. In addition, the values of parameters could be imprecise due to the method of determination employed and a natural variation in values can occur in a population. Consequently, there is some uncertainty regarding the reliability and fidelity of the model predictions and their interpretations. For the purpose of this study we interpret reliability as the extent to which model predictions can be trusted in the context of parameter uncertainty, and fidelity as the degree to which model predictions reflect reality taking heterogeneity in a population into account.
Various methods exist to quantify the effects of parameter variations or changes within a model. These can be used to determine the contribution of parameter uncertainty to variance in model prediction, to test model robustness against parameter changes and even assist in elucidating biologically relevant constituents for intervention with a certain outcome in mind.
In a biological system such as the human body, there is an expected range in which population sizes of cells may fall, as a natural variance will be seen in a population of individuals with slightly different characteristics (i.e. parameter sets). Robustness analysis can demonstrate the possible model outcomes for individuals in a population where biological parameters can vary greatly between individuals. A robust model shows resilience to changes in model inputs, presenting a more stable model [38, 39], although it should be able to account for variances seen in a population. If clinical ranges of observables (such as cell types) are available this analysis can also establish a degree of confidence in the fidelity of model predictions given the variation seen in a population of heterogeneous individuals.
Considering the experimental uncertainty in parameter values, uncertainty analysis allows one to quantify the contribution of uncertainty of a parameter to the overall uncertainty in model predictions [40]. Parameters are often obtained from literature where clinical measurements were pooled in studies not designed for model parametrisation. Uncertainty analysis can thus indicate which parameters have a large effect on model outputs when considering their variances.
Local sensitivity analysis is closely associated to uncertainty analysis, as it entails determining the change in model outputs (e.g. steady state values of model variables) when the inputs or parameters are varied one at a time in a localised parameter space around the reference state [41]. The method is applied to each parameter individually, while the rest are kept at the wild type values, and the results, shown as sensitivity indices, quantify the effect of each parameter on model outcomes near the reference state. Whereas uncertainty analysis indicates which parameters lead to the most uncertainty in model outputs given their variances, local sensitivity analysis quantifies the sensitivities of model outputs to small parameter perturbations around a specific point in parameter space. For a review of possible applications of this analysis see [42, 43]. In the context of metabolic systems, metabolic control analysis (MCA) is a form of local sensitivity analysis which entails calculating normalised partial derivatives of the model outputs (e.g. steady state concentrations or fluxes) with respect to system properties (rates, concentrations or parameters) [44,45,46,47]. In MCA the response coefficient is a sensitivity index quantifying the fractional change in the steady state outputs of the model variables or fluxes upon a 1% change in a parameter [47]. Beyond model sensitivity characterisation, this approach can also be used to indicate possible drug targets and their systemic effects [48,49,50] if one considers parameters with large responses to be potential weaknesses in the system.
The last type of analysis of interest here is a variation on robustness analysis where local sensitivity coefficients are calculated at each point in the complete parameter space. Ideally the parameter ranges should be defined according to the ranges observed in a population. Random parameter sets are constructed and the model simulated for each set  similar to the robustness analysis described above. The local sensitivities (response coefficients) are determined for every parameter set (point in parameter space) and the results are pooled to visualise the spread of response coefficients (sensitivities) in a population. For this analysis the results can be visualised using histograms, and good robustness of local sensitivities in a population is inferred when the results are well conserved with the most probable response coefficient corresponding to the wildtype result from local sensitivity analysis. This analysis can hence assist in determining if the results of the wildtype local sensitivity analysis are approximately retained in a heterogeneous population [48]. As this analysis samples the parameter set from the global parameter space, it could be considered a form of global sensitivity analysis, but since the aim is to analyse the robustness of the local sensitivity results, it will be defined as local sensitivity robustness analysis in this study. Global sensitivity analysis methods used on biological models as described in [47] have different levels of computational and mathematical complexity. It would however be best to use methods that closely relate to the methods used for robustness and local sensitivity analyses for comparison.
This study focusses on analysing four published models of malaria infection where the immune system’s response is incorporated. The models were chosen based on their ability to describe the core dynamics of the disease with varying complexity of the immune system description. Additional considerations were their ordinary differential equations (ODE) structure, and comparability of model variables and processes. Uncertainty, robustness and sensitivity analyses were performed on these withinhost models to determine the effect of parameter uncertainty and variability on the predicted disease dynamics, and to test whether the models could still give reliable and realistic predictions while accommodating heterogeneity and uncertainty.
Model descriptions
The model of Anderson et al. [21] is one of the earliest models on which many others have built. This model includes four variables, describing the RBC, iRBC, merozoite and Tlymphocyte populations, where the Tlymphocyte population represents the immune effectors of the model. The model of Li et al. [22] has a very similar, albeit expanded model structure to that of Anderson et al. [21]. It includes immune effector parameters in Michaelis–Menten–Monod functions to ensure saturation of processes corresponding to the immune system’s response. Niger and Gumel [23] extended the model of Anderson et al. [21] by partitioning the immune system response into two variables: the collective immune effectors and the antibodies specific to malaria infection. The model furthermore separated the iRBC population into different compartments, to account for parasite growth. The final model included for analysis is from the doctoral thesis of Okrinya [24]. Whilst also splitting the immune effectors into two groups as in the Niger and Gumel [23] model, the two immune response variables in this model denote innate and adaptive immunity. It includes an additional variable representing the gametocyte population within the host. The model structures are explained in the following section, and parameter definitions used to reproduce the models prior to analysis can be viewed in Additional file 1: Tables 1–4. It should also be noted that multiple model outputs describing different disease states were obtained in some publications, where parameter values were altered to showcase either parasite free or endemic states. As we are interested in investigating the parameter effects on disease, only models describing infection were used.
Anderson model
Anderson et al. [21] formulated a withinhost model of blood stage malaria infection, which includes the immune response to free roaming merozoites and iRBCs. The model construct with RBCs (x), iRBCs (y), merozoites (s) and Tlymphocytes (T) follows:
The model shows the natural birth rate of healthy RBCs, \(\lambda\), and natural death rates \(\mu x\), \(\alpha y\), ds and aT for RBCs, iRBCs, merozoites and Tlymphocytes respectively. \(\beta {x}{s}\) is a transfer term present in three equations, where \(\beta\) denotes the probability of a merozoite infecting a healthy RBC. Thus, this term depends on, and influences merozoite (s), as well as available RBC population densities (x). The term \(\alpha {r}{y}\) describes the increase of merozoites due to the death rate of iRBC, where the cells burst and release r number of merozoites in the blood. The immune system component decreases the iRBC and merozoite densities with rate constants g and h respectively, whilst also being activated by the iRBCs and merozoites with rate constants k and \(\gamma\). Although this model specifies that Tlymphocytes are the immune effectors used to fight infection, they affect both merozoites and iRBCs. This indicates humoral and cellmediated immunity with humoral immunity fighting against merozoites and cellmediated immunity defending against iRBCs. This model, however, does not include an innate immune response as there is no natural birth rate for immunity effectors. Four different submodels were published based on this model where the last model as shown here uses the complete structure where humoral and cellmediated immunity is included.
Li model
Li et al. [22] contains four equations for variables H (RBCs), I (iRBCs), M (merozoites) and E (immune effectors). Here “immune effectors” include all biological immune effectors and no distinction is made between innate and adaptive immunity.
This model extends the Anderson model, by incorporating nonlinear bounded MichaelisMentenMonod functions to account for saturation of the immunelinked elimination processes. The first process of this kind is the removal of the iRBCs (I) by immune effectors E, represented by \(\frac{{p_1}{I}{E}}{1+\beta {I}}\) in Eq. 2b. In this term \(p_1\) is a rate constant for the rate at which immune effectors can remove iRBCs, and \(\frac{1}{\beta }\) is viewed as a halfsaturation constant for iRBCs. The same holds true for the term \(\frac{{p_2}{M}{E}}{1+\gamma {M}}\), where \(p_2\) is the rate at which the immune effectors can remove the merozoites in the blood plasma, while \(\frac{1}{\gamma }\) is a halfsaturation constant. \(k_1\) and \(k_2\) describe the proliferation rate of lymphocytes due to activation by iRBCs and merozoites, respectively. The immune response is split into two components in the relevant equation 2d. Here the second term corresponds to the immune response component that proliferates due to the activation by iRBCs. The third term corresponds to the immune response component that is activated by merozoites. The merozoites are eliminated by the humoral immune effectors (rate constant \(p_2\)) and the iRBCs by the cellmediated immune effectors (rate constant \(p_1\)). Immune effector activation is a saturable process and depends on the population density of the disease variables as well as the immune cell concentration, their binding and the efficiency of the process. Six submodels were published with the same structure but changes were made in key parameters to obtain different disease dynamics. The one chosen here for analysis describes an endemic state of malaria infection where the immune system’s response is included.
Niger model
Niger et al. [23] contains age compartments (\(Y_i\)) of the intracellular parasite stage (iRBCs), representing the different stages of parasite growth in an infected erythrocyte. The authors also proposed a more physiologically realistic model which splits the immune effectors into two groups: immune cells B and antibodies A. The model, which includes healthy RBCs (X) and merozoites (M), follows:
For the first stage of infection \(Y_1\), Eq. 3b shows infection of healthy RBCs dependent on \(\beta\), the infection rate constant, as well as the variables X and M, depicting the healthy RBCs and merozoite population, respectively. The natural death rate \(\mu _i{Y_i}\) of iRBCs and the progression rate to the next compartment, \(\gamma _{i}{Y_{i}}\), both decrease the \(Y_i\) population in each compartment. The immune system now additionally kills iRBCs, where \(k_i\) is the immunosensitivity of the stage i iRBCs to immune effectors B. In this model n covers five stages where the parameter values differ slightly as shown in Appendix Table 3. B accounts for immune cells including the innate immune effectors, whereas A accounts for merozoites specific antibodies. This distinction is important as innate immune cells are everpresent, whereas antibodies are only formed when the acquired immune response is activated. For the immune cells B there is a natural background production rate of cells, \(\lambda _B\), as well as a stimulation of production rate due to the presence of an infection (at a rate \(\rho _i{Y_i}\)). This stimulation of the production happens due to all infected compartments including free merozoites. For antibodies there is only an increase of antibodies at a rate \(\eta {B}{M}\), dependent on the immune effector and merozoite populations. However, in this model the antibodies (A) are not included in any other processes and is therefore not integrated into the immune response that influence infection. The model used here shows a stable endemic state with an immune response.
Okrinya model
The Okrinya model [24] has an extra variable G for gametocyte population. The model with uninfected RBCs (X), iRBCs (Y), merozoites (M), innate immune cells (P) and adaptive immune cells (A), is given by:
The effector variable A represents the population of adaptive immune cells and include cells like B and Tcells, but the author interprets it as the concentration of antibodies. The immune effectors are thus split differently compared to the Niger model. Comparison to the structure of the Niger model also shows the differences in processes. For example, terms that include \(\frac{1}{(1+c_1{A})}\) indicate the effect of antibodies blocking the infection. The term \({k_m}{P}{M}(1+k_b{A})\) depicts the immune system’s successful removal of merozoites, which is dependent on the concentrations of both the innate immune cells and merozoites. An extra parameter \(\theta\) in the equations for M and G (Eq. 4c and 4d) is the fraction of merozoites that will develop into gametocytes and therefore not reenter the erythrocytic cycle. As was seen for the Niger model, innate immune cells have a natural birth rate and a stimulation rate bought on by infection (\(\eta _1\)), whereas the adaptive immune effectors are only produced when there is an actual presence of disease. Note that a time delay (\(td_1\)) is included in the equation, to account for the interval between start of infection and production of adaptive immune cells. This time delay affects the dynamics upon infection but does not influence the achieved variable steady state values, and was ommited in our reproduction to simplify analysis. The model, describing pathogenesis in an infected individual with an immune response, was published in a dimensional and nondimensionalised form. To enable the comparison of this model to the others, the dimensional model was used.
General remarks on model reproduction
Models were reproduced from literature using the published parameters values as the reference/wildtype set. Simulation results for steady state values of variables and dynamic behaviour were compared to published results and showed good agreement if not exact. It should be noted that for the Li model, some parameter sets used in a Monte Carlo simulation presented here, resulted in oscillatory behaviour (limit cycle). In these cases, our analyses were performed using the unstable steady state.
Parameter symbols and definitions differ between models. Biologically similar parameters were therefore first identified and for the purposes of this study they are denoted by the same symbol (which might differ from the original symbols used in the models) as shown in Table 1. All parameter definitions and values for each model can be viewed in Additional file 1: Tables 1–4.
Results
Robustness analysis
Figure 1 shows the range of endemic steady state values for RBCs and iRBCs for 10,000 parameter sets obtained by a Monte Carlo (MC) random samplingbased technique. All parameters were randomly varied within \(10\%\) of their wildtype values and the steady state values determined for each parameter set. The results for the models of Anderson and Niger show three different disease states i.e. a healthy patient (green), an infected patient with no immune response (red) and an infected patient with an immune response (orange).
Note that the healthy RBCs and iRBCs are unitless in Fig. 1, as Anderson published units only in 1/day and Niger published units in “concentration”. With no infection (green) the healthy RBCs achieve the highest steady state values. These values are shown to drastically decrease with infection (red) and then increase with the addition of the immune response (orange). Although the values of the steady states reached for the infected population with the inclusion of the immune system is higher than those for infection without an immune response, the uninfected steady state values are not recovered. For the iRBCs there are no steady state values in an uninfected individual, nonzero steady state values when there is infection, and subsequently decreased steady state values with the added immune response. These results were shown to indicate the differences between disease states but as is evident, the values reached at steady state for the healthy and infected RBCs in the Anderson model are extremely low. These plots are therefore only to be used as a general representation of how these variable distributions will shift upon infection and the inclusion of the immune response. The general trend of the healthy and infected RBCs as seen for Fig. 1, where, for example, the healthy RBCs decrease with infection and then increase with the inclusion of the immune response, was present in all models with the simulation of the different disease states.
Table 2 shows the steady state value ranges obtained for all variables in each of the models with infection and an immune response. Note that immune effectors include cells from the innate and adaptive immune system. Antibodies (A) from the Niger model is not included in the immune effectors shown here as they are not classified as cells and the immune effectors for the Okrinya model includes only the innate immune cells. Robustness analysis should give all models 10,000 parameter sets with an equivalent number of steady state values per variable, however, some of the Li parameter sets resulted in unphysical steady state values (negative concentrations), possibly due to the instability mentioned in the “General remarks on model reproduction” section above, and were excluded from the results. To compare model outputs using the same number of parameter sets, the Li model therefore required additional sampling to obtain a total of 10,000 valid parameter sets.
Of all the models the Li model showed the best correlation to known clinical ranges for the healthy RBC population in infected individuals (4.41–6.48\(\times 10^6 \;{\mathrm{{cells}}}/\upmu {\mathrm{{L}}}\)) [51]. Immune effectors in the models do not corresponded to the published range of white blood cells (4.52–7.99\(\times 10^3 \;{\mathrm{{cells}}}/\upmu {\mathrm{{L}}}\)) [52]. It is also noteworthy that the range for the Li model’s immune effectors is very large, showing that the immune effectors are very sensitive to changes in the model parameters, which could also affect all other variables. The \(\%\) parasitemia ranges can be calculated from these results using \(100\cdot \frac{iRBCs}{iRBCs+RBCs}\), for simulation results from the same parameter set. Given that severe malaria is classified as \(>5\%\) parasitized cells [53, 54], all models indicate low to moderate parasitemia with the Niger model verging the closest to severe malaria.
For the robustness results in Fig. 2, each model’s results were normalized to their median as all four models achieved different ranges of steady state values.
Although the parameter sets are obtained from uniform distributions, robustness analysis results show more of a normal distribution, indicating good robustness for all models. The Anderson model appears to be the most robust for the merozoites, with the Li model showing a wider range for both the merozoites and iRBCs. The Okrinya and Niger models also show good robustness for the iRBCs and merozoites.
Uncertainty analysis
The uncertainty analysis results indicate which parameters contribute the most to uncertainty in model outputs given their variances. The results are summarized in Table 3 and shows which parameters contribute the most to uncertainty in the merozoite and iRBC populations in all four models.
For the first three models, parameters related to the death rate of immune effectors, \(\mu _P\), and the proliferation rate of immune effectors due to the activation by iRBCs, \(k_1\), contribute the most to uncertainty in merozoites and iRBCs model outputs. These parameters have a larger combined total contribution in the Li model as the total percentage adds to more than \(70\%\), while it is considerably lower in other models. This indicates that there are other parameters in the Anderson and Niger models that can contribute significantly to uncertainty in variable outputs. \(\mu _P\) and \(k_1\) are significant contributors to uncertainty as variance in these parameters can greatly influence the strength of the immune system’s response to infection. In the Okrinya model the number of merozoites released per bursting iRBC and the birth rate of healthy RBCs lead to the most uncertainty in the merozoite population. The more merozoites released per bursting iRBC, the higher the infection and successive infections, leading to large changes in the model outputs. Changes in the birth rate of healthy RBCs (\(\lambda _X\)), lead to changes in the number of RBCs for infection, also leading to significant changes in the model outputs. The death rate of iRBCs (\(\mu _Y\)) similarly has a significant effect due to its impact on parasite production and immune response.
Overall, for three of the four models, the death rate of immune cells and their proliferation rate due to iRBCs have the largest contributions to uncertainty in model outputs. This indicates that the immune system components and lack of detail, contribute the most to uncertainty in the outputs of disease variables.
Local sensitivity analysis
The local sensitivity analysis results concerning comparable parameters between the four models is shown as a heatmap in Fig. 3. Only the disease variables are indicated as it is these variables that are of clinical interest. The results indicate the percentage change in the value of a variable at steady state upon a \(1\%\) change in a parameter value.
Four striking variableparameter pairs are present in the Li model. The death rate of immune effectors, \(\mu _P\), increases both disease variables the most. With a response coefficient larger than 2, the analysis illustrates that a \(1\%\) increase in the death rate of immune effectors will increase the steady state population of merozoites and iRBCs by more than \(2\%\). This is due to a higher death rate of immune cells leading to fewer active immune cells to fight the merozoites. The proliferation rate of immune effectors due to the activation by iRBCs (\(k_1\)) show the exact opposite, where a \(1\%\) increase correlates with a decrease in the disease variable populations by more than \(2\%\). This demonstrates that a small increase in how well the immune system’s response is activated, can lead to a dramatically better disease clearance. Additionally, the response coefficients obtained for the immune effector variables in the Li model showed extremely high responses on parameter perturbations, but was not included in the results shown here as it is a comparison of the disease variables only.
The death rate of immune effectors also delivers high response coefficients for the Anderson and Niger disease variables. These results emphasize the immune system’s response since maximizing the activation of the immune system relative to cell death assists with disease clearance. In the Okrinya model, the number of merozoites released per bursting iRBC, r, has the largest influence on the merozoite population, where the birth rate of healthy RBCs, \(\lambda _X\), influences the iRBCs the most. Logically, the merozoite population will increase if more merozoites are released when a cell bursts and with a larger healthy RBC population, more cells can be infected to form iRBCs.
Interestingly, when inspecting the different models, the Okrinya model shows a direct correlation between the merozoite and iRBC populations, where if a parameter increases the population of merozoites, the iRBC population will also increase. This proportionality in the disease variable responses is also seen with most of the prominent results such as the responses with respect to \(k_1\) and \(\mu _P\) in all models. However, when considering the results for \(\mu _M\) (the natural death rate of merozoites), it is evident that the merozoites decrease whereas the iRBC slightly increase in the models of Anderson and Li. One would expect that if the death rate of merozoites increases and the merozoite population decreases, there would be fewer merozoites to infect RBCs and thus fewer iRBCs. The same should hold for the number of merozoites released per iRBC (r), where an increase in the parameter increases the merozoites and the iRBCs in the Niger and Okrinya models. The counterintuitive results in the Anderson and Li models lay bare the differences in the immune systems in the models, where small changes in infection close to the reference state are responded to differently.
Overall, the local sensitivity analysis revealed the greatest sensitivities for parameters that affect the immune system, the number of merozoites released by bursting iRBCs and the availability of healthy RBCs for infection. Of these processes, the death rate of immune cells and the proliferation rate of the immune effectors due to iRBCs have the largest effects on the disease variables.
Local sensitivity robustness analysis
Using the parameter sets from the robustness analysis, local sensitivity analysis was performed with respect to each parameter, with the model output generated from every parameter set used as a reference state. Response coefficients for a variableparameter pair from all the parameter sets were pooled and visualized using probability distribution histograms. In Fig. 4, the same variableparameter pair is displayed for all models. It illustrates the response coefficient distribution for the response of the iRBCs (Y) to the probability of infection of RBCs by free roaming merozoites. This variableparameter pair was chosen for comparison, as it is also present in the local sensitivity analysis results section.
From the results it is evident that there is an approximately normal distribution for all four models. Furthermore, the response coefficients determined for all of the parameter sets are in a very small range, showing substantial conservation despite heterogeneity in a population setting. Results from the Anderson, Niger and Okrinya models displayed approximately normal distributions for all their variableparameter pairs and their wildtype response coefficients correspond well with the most probable response in a population. The Li results, however, showed some irregularities with the immune effector variable pairs as can be seen in Fig. 5.
It is noteworthy how large the wild type response coefficients are for the Li model with a value of approximately 78, indicating that a \(1\%\) change in the infection rate of healthy RBCs with merozoites, will increase the immune effector population by \(78\%\). The results furthermore reveal the great difference between the wild type, indicated as the black dashed line, and most probable, the peak of the histogram, response coefficients. The wildtype response coefficients obtained with local sensitivity analysis showed large values for the immune effector E responses with respect to changes in most of the parameters. In Fig. 5 it is evident that the most probable response coefficient is much lower given these parameter ranges. Similar local sensitivity robustness behaviour was observed for all parameters that lead to large responses in the immune system variables in the local sensitivity analysis results. It should be noted, however, that the authors investigated the stability of their model and determined bifurcation parameters that produced oscillations. An example of one of the parameters used is \(k_1\), where a value of \(4.5001\times 10^{5}\) resulted in stable steady state values. As this value increases to larger than \(4.5045\times 10^{5}\), a periodic solution appears [22]. These values show how a very small perturbation in some parameters could lead to a large alteration in the model behaviour. Although our analysis always considers the steady state value (or unstable steady state value in the case of oscillations), different parameter sets could lead to greatly varying response coefficients.
Overall, these results suggest that Anderson, Niger and Okrinya models all display conserved sensitivities and responses when considering population heterogeneity, whilst caution should be taken when interpreting analysis results of the Li model.
Discussion
Models describing the immune response to malaria infection were analysed to study the differences in model formalism and sensitivity to parameter values. The methods used in this investigation can indicate when models are more relevant for describing disease dynamics within an individual as well as in a population, whilst also giving insight into the necessity for parameter value certainty.
Robustness analysis was performed to obtain an indication of the spread of outcomes in a population and how it compares to the wild type outcome. Results indicated that the Anderson model was the most robust of all the models showing the most resistance to changes in parameter values and yielding results showing the least variation in outputs. However, this model is the oldest of the models, with limited detail of the interactions between the host immune system and the malaria parasites. The second and third most robust models are those of Okrinya and Niger. Between the two, the Okrinya is more mechanistically descriptive of the withinhost dynamics of malaria infection. Robustness as defined here should, however, be considered in the context of parameter distributions observed in a population since one would expect the model to display a range of predictions that correlate to what is observed in reality given such distributions.
Uncertainty analysis can be helpful to determine which parameters have a large influence on the uncertainty in model outputs. Parameter values were all varied within \(10\%\) around their wild type values for analysis. However, the parameter variances can be much larger in reality, owing to the difficulty of estimating parameter values, as well as variations between individuals due to differences in factors such as diets, ages etc. Gaining more information on these parameters and their possible variances in a population could consequently be useful for further analysis and interpretations. Uncertainty analysis demonstrates how necessary it is to determine precise values or ranges for some parameters. In this study it was found that parameters of the immune system contributed the most to uncertainty in model outputs.
Local sensitivity analysis indicated that the death rate of immune effectors, the number of merozoites released per bursting iRBC and the proliferation of immune cells due to the presence of iRBC have the largest influence on the disease variables in all models around the published reference state. The Okrinya model further indicated the birth rate of healthy RBCs as a parameter of interest. We also found that a larger ratio of proliferation to death rate of the immune cells can drastically affect the model outcomes and help with the attack on infection, and moreover it seems to be the largest influencer in all the models. Combined with decreasing the number of merozoites released per bursting iRBC, this could lead to disease clearance. The two parameters corresponding to the immune system were, however, also present in most of the models as the parameters that contributed the most to uncertainty in model outputs.
Local sensitivity robustness analysis demonstrated how well response coefficients are conserved over the multidimensional parameter space. Good robustness of local sensitivities was achieved for most response coefficients in the models, i.e. the wildtype responses were well conserved in a heterogeneous population. The results of the Li model indicated poor conservation as the large response coefficients observed in local sensitivity analysis did not correspond with the much lower most probable response coefficients. Furthermore, in the Li model the range of some coefficients in a population varied from just above zero to approximately 80, showing the wide range of possible, but intuitively implausible, response coefficients in a population. The results for the Li model thus emphasises the necessity of this analysis, as it demonstrates how local sensitivity analysis may fail in giving insight into the responses in a population of different individuals. As is the case for this model, it could be as a result of transitions in qualitative behaviour such as the switch from steady states, to limit cycle oscillations observed here.
Conclusion
The results from local sensitivity analysis emphasize the importance of the immune response on the disease dynamics of malaria infection, as well as highlight parameters like the death and proliferation rates of immune effectors that could be investigated for disease eradication as these have the largest effect on disease variables. To further establish trust in reliability and fidelity of model predictions the uncertainty in specific parameter values needs to be minimised and clinical ranges obtained where possible. The robustness analysis indicated that the Niger and Okrinya models were robust in the sense that plausible outputs were obtained close to the reference state and local sensitivity robustness analysis showed that the local sensitivity results are relatively well conserved in a population.
In the current study the goal was to establish a measure of confidence in models that encapsulate the core dynamics of a disease state, but the analyses described here can easily be extended to any deterministic models that have similar ODE structures, for malaria and for other diseases. For example, an improved withinhost model could be constructed and analysed with a combination of components from the Niger and Okrinya models where models such as that from Song et al. [55] can be incorporated to include variables of resistant and nonresistant parasites as well as treatment parameters.
Although not the focus of this study, an additional benefit of sensitivity analysis as conducted here is that it can also point to possible drug targets and whether these targets are conserved in a population. Local sensitivity analysis indicated that the number of merozoites released per bursting iRBC has a large effect on all models and that an increased ratio of immune effector proliferation to death rate could be beneficial to disease clearance. The conservation of these results in the local sensitivity robustness analysis shows that this could hold true in a population as well.
Methods of analysis
All reproductions and analyses of the models were completed in Wolfram Mathematica version 12 using standard Mathematica functions. Code used for analyses is provided as executable notebooks in Additional files 2–5, and in pdf format in Additional file 6.
Robustness analysis
Robustness analysis was used to determine the ranges and distributions of the disease variable outputs of each model. This method of analysis can indicate the expected steady state values of withinhost variables in a population and can also show which models are the most resilient to changes in model inputs (i.e. parameters). Monte Carlo (MC) random sampling from uniform distributions was therefore used to obtain 10,000 parameter sets. In the absence of known ranges of parameters and to facilitate direct model output comparisons, all parameters were varied simultaneously within a \(10\%\) range of their wild type values for all models. These sets mimic different parameter sets of individuals within a population. All sets were used to determine the distribution of the disease variables’ values at steady state. Three different disease states were simulated to determine the difference in the ranges of the steady state values for the RBC and iRBC populations: (i) no infection, (ii) infection without an immune response and, (iii) infection with an immune response. Results for two of the models are shown in the results section for comparison. The submodels where the immune system is incorporated from each of the four publications were then analysed with regards to their disease variables in the endemic state. Where required for comparison between models, in silico datasets were normalized to their median. The results are visualized with boxandwhisker plots.
Uncertainty analysis
Here we followed the method set out in [56]. To determine the uncertainty in variable outputs due to variability in parameter values, the variance (\(\sigma ^2\)) of the natural logarithm of each parameter (\(p_j\)) is first calculated:
However, as the variance of many of the parameters in these biological models are not known, the variance of each parameter was calculated using a uniform distribution around its wildtype value with upper and lower bounds given by wildtype value \(\pm 10\%\). The results therefore indicate which parameters contributes the most to uncertainty in the model variable outputs. The individual contribution of each parameter variance \(\sigma ^2_j(\ln {p_j})\) to the total variance of each variable \(\sigma _j^2 (V_i )\) is then calculated using:
The total variance of each variable can be determined by summation of the individual contributions:
The contribution of each parameter to model variable uncertainty is then given by:
Local sensitivity analysis
To determine for which parameters the disease variables are most sensitive in the reference state, local sensitivity analysis was performed . The method entails the perturbation of one parameter at a time to see the effect on model outputs, indicated as a response coefficient. A response coefficient describes the percentage change of a model output — in this case the steady state values of different variables V — upon a \(1\%\) change in model inputs or parameters p:
Numerically, the derivative is approximated by a finite difference formula using small perturbations around the wild type parameter value and noting the change on model output.
Local sensitivity robustness analysis
Local sensitivity robustness analysis was used to test for conservation of the local sensitivity analysis results given parameter variations in a population. It would therefore be an indicator of whether the response coefficient results of a parametervariable pair in the individual with the wild type parameter set, is the most common result in a population with differing parameter sets. Local sensitivity analysis, as described above, was therefore performed for 10,000 parameter sets obtained using MC random sampling from the uniform distributions of parameters in the complete parameter space (see “Robustness analysis” above). The response coefficients for a given parameter was pooled from each set and histograms were used to visualize results.
Availability of data and materials
All models analysed in this study are available in the curated database of JWS Online: https://jjj.bio.vu.nl/models/anderson1, https://jjj.bio.vu.nl/models/niger1, https://jjj.bio.vu.nl/models/li2, https://jjj.bio.vu.nl/models/okrinya1 Model parameter values and variable definitions are also provided in Additional file 1. Code used for analyses is provided as executable notebooks in Additional files 2–5, and in pdf format in Additional file 6.
Abbreviations
 RBC:

Red blood cell
 iRBC:

Infected red blood cell
 MCA:

Metabolic control analysis
 MC:

Monte Carlo
References
 1.
World Malaria Report 2019. World Health Organization. https://www.who.int/publications/i/item/9789241565721. Accessed 1 Mar 2021.
 2.
Hanboonkunupakarn B, White NJ. The threat of antimalarial drug resistance. Trop Dis Travel Med Vaccines. 2016;2(1):1–5.
 3.
Cox FEG. History of the discovery of the malaria parasites and their vectors. Parasit Vectors. 2010;3(1):1–9.
 4.
Cowman AF, Healer J, Marapana D, Marsh K. Malaria: biology and disease. Cell. 2016;167(3):610–24.
 5.
White NJ. Review series: antimalarial drug resistance. J Clin Investig. 2004;113(8):1084–92.
 6.
Siciliano G, Alano P. Enlightening the malaria parasite life cycle: bioluminescent plasmodium in fundamental and applied research. Front Microbiol. 2015;6(MAY):1–8.
 7.
Hato T, Dagher PC. How the innate immune system senses trouble and causes trouble. Clin J Am Soc Nephrol. 2015;10(8):1459–69.
 8.
Abbas A, Litchman A, Pillai S. Basic immunology: functions and disorders of the immune system. 5th ed. Amsterdam: Elsevier; 2015.
 9.
Bonilla FA, Oettgen HC. Adaptive immunity. J Allergy Clin Immunol. 2010;125(2 SUPPL. 2):S33–40.
 10.
Deroost K, Pham TT, Opdenakker G, Van den Steen PE. The immunological balance between host and parasite in malaria. FEMS Microbiol Rev. 2016;40(2):208–57.
 11.
SherrardSmith E, Skarp JE, Beale AD, Fornadel C, Norris LC, Moore SJ, et al. Mosquito feeding behavior and how it influences residual malaria transmission across Africa. Proc Natl Acad Sci U S A. 2019;116(30):15086–96.
 12.
Bongfen SE, Torgler R, Romero JF, Renia L, Corradin G. Plasmodium berghei—infected primary hepatocytes process and present the circumsporozoite protein to specific CD8 + T cells in vitro. J Immunol. 2007;178(11):7054–63.
 13.
Murphy SC, Shott JP, Parikh S, Etter P, Prescott WR, Stewart VA. Review article: malaria diagnostics in clinical trials. Am J Trop Med Hyg. 2013;89(5):824–39.
 14.
Kwiatkowski D, Nowak M. Periodic and chaotic hostparasite interactions in human malaria. Proc Natl Acad Sci U S Am. 1991;88(12):5111–3.
 15.
Tabo Z, Luboobi LS, Ssebuliba J. Mathematical modelling of the inhost dynamics of malaria and the effects of treatment. J Math Comput Sci. 2017;17(1):1–21.
 16.
Orwa TO, Mbogo RW, Luboobi LS. Mathematical model for hepatocyticerythrocytic dynamics of malaria. Int J Math Math Sci. 2018;2018:1–18.
 17.
Dietz K, Raddatz G, Molineaux L. Mathematical model of the first wave of Plasmodium falciparum asexual parasitemia in nonimmune and vaccinated individuals. Am J Trop Med Hyg. 2006;75(SUPPL. 2):46–55.
 18.
Nannyonga B, Mwanga GG, Haario H, Mbalawata IS, Heilio M. Determining parameter distribution in withinhost severe P. falciparum malaria. BioSystems. 2014;126:76–84.
 19.
Diebner HH, Eichner M, Molineaux L, Collins WE, Jeffery GM, Dietz K. Modelling the transition of asexual blood stages of Plasmodium falciparum to gametocytes. J Theor Biol. 2000;202(2):113–27.
 20.
Eichner M, Diebner HH, Molineaux L, Collins WE, Jeffery GM, Dietz K. Genesis, sequestration and survival of Plasmodium falciparum gametocytes: parameter estimates from fitting a model to malariatherapy data. Trans R Soc Trop Med Hyg. 2001;95(5):497–501.
 21.
Anderson RM, May RM, Gupta S. Nonlinear phenomena in hostparasite interactions. Parasitology. 1989;99(S1):S59–79.
 22.
Li Y, Ruan S, Xiao D. The withinhost dynamics of malaria infection with immune response. Math Biosci Eng. 2011;8(4):999–1018.
 23.
Niger AM, Gumel AB. Immune response and imperfect vaccine in malaria dynamics. Math Popul Stud. 2011;18(2):55–86.
 24.
Okrinya AB. Mathematical modelling of malaria transmission and pathogenesis. Loughborough: Loughborough University; 2014.
 25.
Eckhoff PA. Malaria parasite diversity and transmission intensity affect development of parasitological immunity in a mathematical model. Malaria J. 2012;11:1–14.
 26.
McKenzie FE, Bossert WH. The dynamics of Plasmodium falciparum bloodstage infection. J Theor Biol. 1997;188(1):127–40.
 27.
Yan Y, Adam B, Galinski M, Kissinger JC, Moreno A. Mathematical model of susceptibility, resistance, and resilience in the withinhost dynamics between a Plasmodium parasite and the immune system. Math Biosci. 2015;270(Pt B):213–23.
 28.
Molineaux L, Diebner HH, Eichner M, Collins WE, Jeffery GM, Dietz K. Plasmodium falciparum parasitaemia described by a new mathematical model. Parasitology. 2001;122(Pt 4):379–91.
 29.
Hellreigel B. Modelling the immune response to malaria with ecological concepts: shortterm behaviour against longterm equilibrium. Proc Biol Sci. 1992;250(1329):249–56.
 30.
Tumwiine J, Mugisha JYT, Luboobi LS. On global stability of the intrahost dynamics of malaria and the immune system. J Math Anal Appl. 2008;341(2):855–69.
 31.
Chiyaka C, Garira W, Dube S. Modelling immune response and drug therapy in human malaria infection. Comput Math Methods Med. 2008;9(2):143–63.
 32.
Anderson RM. Complex dynamic behaviours in the interaction between parasite populations and the hosts immune system. Int J Parasitol. 1998;28(4):551–66.
 33.
Gurarie D, McKenzie FE. Dynamics of immune response and drug resistance in malaria infection. Malaria J. 2006;5:1–15.
 34.
Klein EY, Graham AL, Llinás M, Levin S. Crossreactive immune responses as primary drivers of malaria chronicity. Infect Immun. 2014;82(1):140–51.
 35.
De Leenheer P, Pilyugin SS. Immune response to a malaria infection: properties of a mathematical model. J Biol Dyn. 2008;2(2):102–20.
 36.
Gurarie D, Karl S, Zimmerman PA, King CH, Pierre TG, Davis TME. Mathematical modeling of malaria infection with innate and adaptive immunity in individuals and agentbased communities. PLoS ONE. 2012;7(3):1–13.
 37.
Antia R, Nowak MA, Anderson RM. Antigenic variation and the withinhost dynamics of parasites. Proc Natl Acad Sci U S A. 1996;93(3):985–9.
 38.
Zhu Y, Wang QA, Li W, Cai X. Analytic uncertainty and sensitivity analysis of models with input correlations. Phys A Stat Mech Appl. 2018;2018(494):140–62.
 39.
Rizk A, Batt G, Fages F, Soliman S. A general computational method for robustness analysis with applications to synthetic gene networks. Bioinformatics. 2009;25(12):169–78.
 40.
Schuster S, Heinrich R. The definitions of metabolic control analysis revisited. BioSystems. 1992;27(1):1–15.
 41.
Charzyńska A, Nalńcz A, Rybiński M, Gambin A. Sensitivity analysis of mathematical models of signaling pathways. Biotechnologia. 2012;93(3):291–308.
 42.
Zi Z. Sensitivity analysis approaches applied to systems biology models. IET Syst Biol. 2011;5(6):336–46.
 43.
Borgonovo E, Plischke E. Sensitivity analysis: a review of recent advances. Eur J Oper Res. 2016;248(3):869–87.
 44.
Heinrich R, Rapoport TA. A linear steadystate treatment of enzymatic chains. General properties, control and effector strenght. Eur J Biochem. 1974;42(1):89–95.
 45.
Hofmeyr JHS. Metabolic control analysis in a nutshell. In: Hucka M, Morohashi M, Kitano H, et al, editors. Proceedings of the 2nd International Conference on Systems Biology; 2001 Nov 57; California Institute of Technology, Pasadena, USA. Wisconsin: Omnipress; 2001. p. 291–300.
 46.
Kacser H, Burns JA. The control of flux. Symp Soc Exp Biol. 1973;27:65–104.
 47.
De Pauw DJW, Vanrolleghem PA. Practical aspects of sensitivity function approximation for dynamic models. Math Comput Model Dyn Syst. 2006;12(5):395–414.
 48.
Jurina T, Tušek AJ, Čurlin M. Local sensitivity analysis and metabolic control analysis of the biological part of the BTEX bioremediation model. Biotechnol Bioprocess Eng. 2015;20:1071–87.
 49.
Saitou T, Itano K, Hoshino D, Koshikawa N, Seiki M, Ichikawa K, et al. Control and inhibition analysis of complex formation processes. Theor Biol Med Model. 2012;9(1):1.
 50.
Özbayraktar FBK, Ülgen KÖ. Drug target identification in sphingolipid metabolism by computational systems biology tools: metabolic control analysis and metabolic pathway analysis. J Biomed Inform. 2010;43(4):537–49.
 51.
Omuse G, Maina D, Mwangi J, Wambua C, Radia K, Kanyua A, et al. Complete blood count reference intervals from a healthy adult urban population in Kenya. PLoS ONE. 2018;13(6):1–19.
 52.
Kotepui M, Phunphuech B, Phiwklam N, Chupeerach C, Duangmano S. Effect of malarial infection on haematological parameters in population near ThailandMyanmar border. Malaria J. 2014;13(1):1–7.
 53.
Tekeste Z, Workineh M, Petros B. Determining the severity of Plasmodium falciparum malaria in Ethiopia. J Infect Public Health. 2013;6(1):10–5.
 54.
Trampuz A, Jereb M, Muzlovic I, Prabhu RM. Clinical review: severe malaria. Critical Care. 2003;7(4):315–23.
 55.
Song T, Wang C, Tian B. Mathematical models for withinhost competition of malaria parasites. Math Biosci Eng. 2019;16(6):6623–53.
 56.
Zádor J, Zsély IG, Turányi T. Local and global uncertainty analysis of complex chemical kinetic systems. Reliab Eng Syst Saf. 2006;91(10–11):1232–40.
Acknowledgements
Not applicable.
Funding
We acknowledge the financial assistance from the DST/NRF in South Africa, particularly for funding the SARChI initiative (Jacky Snoep NRFSARCHI82813), and the CSUR (David van Niekerk ID 116298) and Innovation Doctoral Scholarships (Shade Horn ID 121296) programmes.
Author information
Affiliations
Contributions
D.D.v.N. and J.L.S. designed research; S.H. performed research; D.D.v.N. and J.L.S. provided resources and funding; D.D.v.N. and S.H. wrote the paper. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
No ethics approval was required for this study.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file
Additional file 1.
Model parameters and variables.
Additional file 2.
Anderson analysis (Mathematica notebook).
Additional file 3.
Li analysis (Mathematica notebook).
Additional file 4.
Niger analysis (Mathematica notebook).
Additional file 5.
Okrinya analysis (Mathematica notebook).
Additional file 6.
Analysis code for all models.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Horn, S., Snoep, J.L. & van Niekerk , D.D. Uncovering the effects of heterogeneity and parameter sensitivity on withinhost dynamics of disease: malaria as a case study. BMC Bioinformatics 22, 384 (2021). https://doi.org/10.1186/s1285902104289z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1285902104289z
Keywords
 Withinhost
 Sensitivity analysis
 Malaria
 Modelling