Sensitivity analysis for a within-human-host immuno-pathogenesis dynamics of Plasmodium falciparum parasites

Sensitivity analysis has become increasingly useful in many fields of engineering and sciences. Researchers use sensitivity and uncertainty analysis in the mathematical modelling of biological phenomena because of its value in identifying essential parameters for model’s output. Moreover, it can help in the process of experimental analysis, model order reduction, parameter estimation, decision making or development of recommendations for decision makers. Here, we demonstrate the use of local sensitivity analysis to understand the influence of different parameters on a threshold parameter, RI 0, resulting from the analysis of a within human-host model for the dynamics of malaria parasites. Our results reveal that the obtained RI 0 is most sensitive to the infection rate of healthy red blood cells (RBCs) by merozoites, the average number of merozoites released per bursting parasitized RBCs, the proportion of parasitized RBCs that continue asexual reproduction and the per capita natural death rate of merozoites. Woldegebriel A. Woldegerima African Institute for the Mathematical Sciences (AIMS) Cameroon P.O. Box 608, Limbe, Cameroon Partly developed while visiting Lehigh University as a Pre-doctoral visiting Scholar Gideon A. Ngwa Department of Mathematics, University of Buea, P.O. Box 63, Buea, Cameroon e-mail: akumhed@yahoo.com Miranda I. Teboh-Ewungkem Department of Mathematics, Lehigh University, Bethlehem, PA, 18015, USA e-mail: ewungkems@gmail.com;mit703@lehigh.edu Citation: Woldegebriel A. Woldegerima, Gideon A. Ngwa, Miranda I. Teboh-Ewungkem, Sensitivity analysis for a within-human-host immuno-pathogenesis dynamics of Plasmodium falciparum parasites, in R. Anguelov, M. Lachowicz (Editors), Mathematical Methods and Models in Biosciences, Biomath Forum, Sofia, 2018, pp. 140-168, http://dx.doi.org/10.11145/texts.2018.05.257 Copyright: c © 2018 Woldegerima et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.


Introduction and background
Malaria remains one of the most prevalent human diseases worldwide and still causes a significant problem in many tropical areas, especially in the tropical African region.The malaria control problem comes with a range of challenges that must be surmounted.These challenges, including the issues of morbidity and mortality, are listed in the world malaria report [32].Mathematical and statistical sensitivity analysis can be used to quantify the magnitude and relative importance of some of the malaria disease parameters towards the malaria control problem.Such a sensitivity analysis will be useful if appropriate models are developed that can simulate the malaria control problem and that can complement theoretical and experimental work.
Sensitivity refers to the degree to which an input parameter influences a model's output, and thus sensitive parameters are those which have a significant influence on the model's outcomes [12].Sensitivity analysis is the systematic investigation of a function, a model or a system's response/output to either perturbation of the quantitative factors (input and/or parameters) or variations in the qualitative factors (structure, connectivity modules or submodels) [16].It can be used to study how the uncertainty in the output of a mathematical system is influenced by different sources of uncertainty in its inputs such as the parameters [22] and in identifying influential model parameters that can aid in optimizing model structure.It can play a crucial role for decision making policy in developing public health strategies for control and prevention of infectious diseases [19,35].The importance and usage of sensitivity analysis in decision making is well documented [4,12,18,19,22,24].
There are generally two classes of sensitivity analysis approaches: Local Sensitivity Analysis (LSA) and Global Sensitivity Analysis (GSA) [8,15,16,29].Local sensitivity analysis, as defined in Zhou et al. [36], is the assessment of the local impact of an input factors' variation on a model's response whereby concentration on the sensitivity is in the vicinity of a set of factor values.On the other hand, global sensitivity analysis is the process of apportioning the uncertainty in outputs to the uncertainty in each input factor over their entire range of interest [37].To conduct a sensitivity analysis, the input and output variables of the model under study must be clearly understood.For example, in disease dynamics, possible useful outputs are the state solutions, the threshold parameters such as the basic reproduction or the effective reproduction numbers, defined using the system parameters.In [15], the author defined a variance-based global sensitivity analysis which focusses on the variance of the model output.More precisely, he analysed how input variability influences output variance, giving a rigorous quantitative overview.We now present a summary of the local versus global sensitivity analyses.Effect of an input/parameter on an output is obtained by varying inputs about a fixed point by a small amount one at a time.
Effects of inputs on an output are obtained by varying all inputs simultaneously over their reasonable parameter ranges [6,21].
Focus is on estimating local impact of a parameter on model output.
Focus is on analysing global impact on the entire parameter space [31].
Provides a particular view on the influence of an individual input/parameter on a model output, holding other parameters fixed.So changes to input are limited to one parameter at a time and so cannot account for the effects on the output of the likely interactions between parameters which can occur in disease models [4,16,38].
Provides an overall view on the influence of inputs on outputs varying all or some parameters simultaneously.Allows for a possibility that all input parameters of the model can be varied simultaneously and the effects on the output of both individual inputs and interactions between inputs are assessed.
Less accurate but much faster and easy to handle.
More accurate, needs more effort, intensive computational procedures and is challenging in general.[37].
Investigates the behaviour of a model/output of a model in a specific set of nominal inputs/nominal parameter values.
Can better handle input values that are uncertain and cover large ranges.
Uses mainly global statistical methods such as Sampling-based methods using Monte-Carlo techniques, Variance-based methods such as the Method of Sobol [16,29,31,35].

Approaches to sensitivity analysis
There are various methods of local and global sensitivity analysis.See, [1,3,4,5,6,8,9,15,16,29,19,31,35], for example.Local sensitivity methods can be broadly divided into Forward Sensitivity Analysis, which involves approximating derivatives, and Adjoint Sensitivity Analysis which involves evaluating integrals.Our focus will be on the Forward method that involves estimating derivatives.In effect, this local method invokes finding the derivative of the output with respect to an input factor, where the derivative is taken at some fixed point in the space of the input.This derivative is called sensitivity coefficient of the output to the input.If the local sensitivity analysis is for an output with a single input, derivatives will be ordinary, and the sensitivity coefficient can be approximated as we now demonstrate.Finite Difference Method: Consider an output given by a real valued function, φ (p), where p is an input parameter.To study how a small perturbation in p, say ∆ p, influences the response in the output, we use the derivative of φ with respect to p taken at the input p.It is often assumed that Φ has a continuous dependence on the parameter p so that we can, for example, expand φ (p + ∆ p) in a Taylor series in p to approximate the first derivative, after neglecting higher order terms in ∆ p, as This is a forward-difference estimate and is a first-order approximation for the first derivative since the truncation error is O(∆ p).On the other hand, the centraldifference estimate for the first derivative is a second-order approximation since the truncation error is O(∆ p 2 ), and is given as Approximations for second-order and higher-order derivatives can be obtained similarly.For instance, computing the central difference approximation Φ (p + 1 2 ∆ p) and φ (p − 1 2 ∆ p), with ∆ p replaced by 1 2 ∆ p, we have Now, in general, physical problems usually involve the interaction of more than one input and are thus modelled by a system of equations.Therefore, in the above discussion, the output function will be a function of several input parameters, and hence to look for local sensitivities we need to use the gradient of the output perturbed about the corresponding components of the vector of parameters.Consider a parameterized model system arising from an applied problem (say for example in ecology or biology etc) with n− state variables and m− input parameters: where is the vector with components describing the right hand side functions of the system.Particularly, for each i = 1, 2, • • • , n, we have the single model equation Next, we use partial derivatives to determine local sensitivities by computing the changes in the model output (i.e.model solutions to the variations in the model inputs (parameters) at local points in the parameter space.The challenge with using finite-difference formulas to estimate sensitivity coefficients is in the step size selection of the input parameter, ∆ p. Questions relating to whether the step size is too small or too big need to be addressed [3,5].
Calculating sensitivities with finite-differences require at least m + 1 runs of the model since it must be calculated for each of the m parameters for x i [1].Addi-tionally, the finite difference technique is often used without consideration of the non-linearity of the model or model outputs) nor the round-off error introduced by the output calculations.Symbolic (Direct) Differentiation Method: More sophisticated methods involve using the symbolic (direct) differentiation method [5,6,10,36].Consider the singlestate equation in (3).Then differentiating both sides with respect to p j , we have Using the chain rule for derivatives and Clairaut's theorem for mixed partial derivatives, we obtain or, equivalently, is the timederivative of the local sensitivity coefficients (LSE), s i j .We can find the local sensitivity matrix, S j (t), for each j = 1, 2 If we apply chain rule on the left hand side and Clairaut's theorem for mixed partial derivatives on the right hand side, we get where, for each Therefore, for each j = 1, 2, • • • , m the local sensitivity matrix, S j , can be obtained from equation (4) at any time t, where J, F j , S j are as given in (5).Thus, to find S j , we need to solve the coupled system defined in equations ( 2) and (4) given as: Generally, taking the derivative of the components of the vectors in system (2) with respect to the components of the vector parameter p, we arrive at We notice here that the initial conditions S j (t 0 ), j = 1, 2, • • • , m for system (4) or S(t 0 ) can be obtained from the initial conditions of X, since S j (t 0 ) = ∂ X ∂ p j (t 0 ).That is, the initial condition of the local sensitivity matrix S j (t) are simply the partial derivatives with respect to the input parameter p j of the initial condition of X.So, if x i (t 0 ) does not contain p j , for at least one i, then s i j (t 0 ) = 0, ∀(i, j).On the other hand, if the initial condition for x l , 1 ≤ l ≤ n contains p j , then s i j (t 0 ) = 1, only when i = l, otherwise s i j (t 0 ) = 0, for all i = l.Most often, we take s i j (t 0 ) = 0, ∀1 ≤ i ≤ n, 1 ≤ j ≤ m, which means S(t 0 ) = 0 [5,29,38].
We have determined the expressions for the local sensitivity coefficients, s i j (t).The question is what do they measure?Now, since, s i j (t) = ∂ x i ∂ p j , it implies that a one unit change in the input p j will produce a change of s i j in the output x i [29].However, since the size of x i will typically impact the size of s i j , we will consider relative changes (in percents).For this, we introduce normalized local sensitivity coefficients (also referred to by some authors as elasticity index or local sensitivity index), see [29,38] defined here as: The normalized local sensitivity coefficient si j is dimensionless and quantifies the relative changes, giving the percent change in the output (which is x i here) due to a 1% change in the parameter (p j here).The direct differential method presented above has it's advantages and drawbacks when compared to the finite-difference method.Comparatively and as an advantage, the difficulty in step size selection is not a main issue in the differential method but tends to cause problems in the finite-difference method [38].Another advantage of the differential method over the finite difference method is that the sensitivity of different outputs (state variables) with respect to certain parameters can be solved simultaneously.This can be done by simultaneously solving the coupled system (6) for both the state variables X and S j = ∂ X ∂ p j .However, from a disadvantageous standpoint, solving the coupled system ( 6) for all j = 1, 2, • • • , m simultaneously is computationally-demanding since it requires that we solve (m + 1)n ODEs m times [29].Furthermore, determining the Jacobian matrix J is time-consuming especially for large-scale problems.Additionally, when several parameters are changed simultaneously, investigating their effects on model results via the sensitivity matrix S, using the direct differential method is also time-consuming and computationallydemanding.To reduce the time-consuming computation cost, an alternative method called adjoint method (some authors called it the Green function method) was developed [11,29].The adjoint method works by defining a Greens or kernel function matrix.For more details on this method, see [5,6,11,20,29,38] and references therein.

Normalized sensitivity analysis applied to a within-human-host model for malaria parasitemia
Here, we investigate the local sensitivity of the basic reproduction number to key associated parameters for a within human-host immuno-pathogenesis model developed and analysed in [34].For the model developed in [34], the basic reproduction number is a parameter-dependent output of the model system and gives an indication of the severity of Plasmodium parasitemia in an individual, and lowering the number below unity or in some cases below a critical number that is less than 1 [25], is the major way to reduce parasitemia and relieve a patient from illness.Thus studying the monotonicity relations between the reproduction number and the parameters used in the model is important.

The model
The model developed and analysed in [34] incorporates seven main interacting actors for the malaria parasite dynamics within a human-host at time t.These are the densities of (i) healthy/uninfected red blood cells (HRBCs), R h (t), (ii) parasitized red blood cells (IRBCs) (these are healthy red blood cells infected by merozoites), R p (t), (iii) free floating merozoites, M(t), (iv) early/immature stage gametocytes (stages I-IV), G e (t), (v) the late/matured gametocytes (stage V), (G l (t)), (vi) the human innate immune system effectors, E i (t), and (vii) the human adaptive immune system effectors, E a (t).The state variables are presented in Table 2 and relate to each other through the equations: In system ( 8) -( 14), HRBCs (R h ) are recruited following a logistic model and some are parasitized due to contact (assumed to be mass action contact) with free floating merozoites (M), producing IRBCs (R p ) but also lead to elimination of the free floating merozoites.Some of the free floating merozoites also seek to infect parasitized (infected) red blood cells (IRBCs).IRBCs follow one of three possible paths -they either die naturally, or burst to release r merozoites per IRBC, or continue the path towards the formation of the early state gametocytes (G e ) which then leads to the transmissible forms of the malaria parasites, the late state gametocytes (G l ).
In some of the processes, the innate (E i ) and adaptive (E a ) immune system cells are activated/excited and seek to control the mechanisms of parasitization, bursting, contacts, and early state gametocyte formation.The different parameters describing transformations and contact rates are described in Table 3, with ω = Λ − µ h the constant growth rate of healthy red blood cells and K = Λ − µ h μh = ω μh > 0, the maximum carrying capacity of healthy red blood cells.For more details on the model derivation, see [34].Mass action contact rate between free merozoites and HRBCs.Models the effective parasitization rate of HRBCs by merozoites.
Adjusted mass action contact rate between free merozoites and HRBCs.Models the effective absorption rate of free merozoites by HRBCs as the merozoites seek to parasitize HRBCs.In the process, free merozoites are cleared from the blood.
Mass action contact rate between free merozoites and IRBCs.Models the effective absorption rate of free merozoites by IRBCs as the merozoites seek to infect red blood cells (in this case IRBCs) In the process, free merozoites are cleared from the blood.
Per capita natural death rates of healthy red blood cells.
Additional death of healthy red blood cells due to density dependent contact inhibition.

Λ
Linear growth rate of red blood cells due to per capita production of red blood cells from the bone marrow.
T −1 µ p Per capita natural linear death rate of parasitized erythrocytes.
Per capita natural linear death rate of immature gametocytes.
Per capita natural linear death rate of mature gametocytes.
Per capita natural linear death rate of freely-floating merozoites.
Per capita natural linear death rate of adaptive immune system cells.
Effective linear growth rate of innate immune system cells.
Effective carrying capacity of the environment for innate immune system cells.

I
ξ 0 Efficiency of adaptive immune effectors in blocking invasion of healthy red blood cells by merozoite during mass action contact.
Efficiency of adaptive immune effectors in inhibiting merozoite transformation in parasitized red blood cells.
Efficiency of adaptive immune effectors in inhibiting maturation of early state gametocytes.
Elimination/killing rate of parasitized red blood cells by innate immune system cells due to mass action contact.
Elimination/killing rate of free merozoites by innate immune system cells due to mass action contact.
Elimination/killing rate of immature gametocytes by innate immune system cells due to mass action contact.
Elimination/killing rate of mature gametocytes by innate immune system cells due to mass action contact.
Mass action contact rate between infected red blood cells, innate immune system cells and adaptive immune system cells leading killing or elimination of IRBCs.These are additional clearances due to the presence of adaptive immunity which enhances the innate.
ρ n Mass action contact rate between merozoites, innate immune system cells and adaptive immune system cells killing or elimination of merozoites.These are additional clearances due to the presence of adaptive immunity.
Antibodies not only blocks invasion but also assist innate cells in engulfing merozoites.
Mass action contact rate between immature gametocytes, innate immune system cells and adaptive immune system cells.These are additional clearances due to presence of adaptive immunity.
Average number of merozoites released by each bursting IRBC.
Average number of early stage gametocytes arising from one IRBC.
is the proportion of the infected red blood cells that differentiate following the path towards gametocytogenesis.

γ p
Rate at which parasitized red blood cells mature to a point where they either burst to release more free-floating merozoites or continue differentiating towards the gametocytogenesis path.
T −1 γ l Rate at which early stage/immature gametocytes mature to become late stage/mature gametocytes.
Linear response/production rate of innate immune effectors due to stimulation by parasitized red blood cells.
Linear response/production rate of adaptive immune effectors due to stimulation by free floating-merozoites.
Linear response/production rate of adaptive immune system cells due to stimulation by parasitized red blood cells.
Linear response/production rate of adaptive immune system cells due to simulation by free-floating merozoites.
Mass action contact parameter between infected red blood cells and innate immune system cells modelling the rate of depletion of the innate immune.
Mass action contact rate between free-floating merozoites and innate immune system cells, modelling the rate of depletion of innate immune cells due to such contact.
Mass action contact parameter between parasitized red blood cells and adaptive immune system cells, modelling the rate of depletion of adaptive immune cells due to such contact.
Mass action contact parameter between free-floating merozoites and adaptive immune system cells, modelling the rate of depletion of adaptive immune cells due to such contact.
The basic properties of boundedness, positivity well-posedness of solutions for the model have been discussed in [34].Here our main objective is to study sensitivity of the basic reproduction number, an output of the system ( 8)-( 14) with respect to the important parameters of the system.We start by computing the expression for the basic reproduction number using the next generation matrix.For this, we need to first determine the infection-free steady state, which for our within-human-host dynamics, is referred to as the Merozoite Free Steady State (MFSS).The steady state solutions, , of system ( 8)-( 14) are solutions of the algebraic system obtained by setting the right hand side to zero.Starting with (8), Substitute equation ( 15) into the next equation where dR p dt = 0.For the case But, the latter is impossible since γ p , µ p , ρ p > 0, ρ a > 0, and 8), ( 9) or ( 10) each set to zero, yield R * p = 0. Thus, M * = 0 if and only if R * p = 0.By a similar argument applied to Additionally, for equations 13) and E * a = 0 from equation (14).Thus, in the absence of the malaria parasites (i.e.absence of disease), equivalent to R * p = 0 = M * , the possible merozoite-free, and in general parasite-free, steady states occur when either (i) Thus we have the following steady states when there are no parasites within human blood (i.e. when M * = 0): Henceforth, we will refer to X (R h =0,E i =0,E a =0) as the trivial steady state (TSS), ) as the merozoite-free steady state (MFSS), and ) as the boundary steady states.
Biologically, the merozoite-free steady state (MFSS), X (R h =K,E i =K i ,E a =0) = (K, 0, 0, 0, 0, K i , 0), depicts a healthy individual who has never been infected with the malaria parasite or was infected but may have been away from a malaria en-demic region for a while and/or has not received a successful sporozoite (transmissible forms of the parasite from mosquitoes to humans [25,26,27,28]) inoculation from an infective mosquito.Thus, this individual has only their healthy red blood cells and innate immune cells.Their adaptive immune cell population is at the zero level because adaptive immune cells are developed and maintained with continuous re-exposure to the malaria parasite, waning with time when there is no contact between the individual and an infectious mosquito.The first boundary steady state (0, 0, 0, 0, K i , 0), depicts a scenario in an individual who has innate immune cells but their healthy red blood cells have been depleted.In other words, such an individual has only white blood cells but no red blood cells.We understand that such a scenario may be rare but can hypothetically occur if there is heavy parasitemia that compromises the individual's system, leading to the depletion of their healthy red blood cells and likely death of the individual, especially if there is no recourse to remedy the situation.The second boundary steady state (K, 0, 0, 0, 0, 0, 0), depicts a scenario in an individual who has healthy red blood cells but no innate immune cells.In other words, the individual has only red blood cells but no white blood cells.This will be an individual with a very compromised immune system.Definition 1.As defined in [14], by the within-human-host basic reproduction number of the malaria parasite, hereafter denoted by R I 0 > 1, we mean the average number of second generation infected red blood cells generated by one primary infected red blood cell at the onset of infection.Thus, if R I 0 > 1, then on average each IRBC produces more than one new IRBC and hence the merozoites are able to invade and infect available HRBCs implying persistence of parasitemia.On the other hand, if R I 0 < 1, then on average an IRBC produces less than one new IRBC and hence the parasite is not able to establish in the human, leading to possible parasite clearance in the human.R I 0 = 1 is a threshold value below which the generation of secondary cases is insufficient to maintain the infection within the human-host.Now, the creation of newly infected red blood cells comes as a result of the bursting of the IRBCs to release of merozoites.Thus, R I 0 as defined above can be thought of as a measure of the efficiency of a merozoite's first infection.
Next, we use the next-generation matrix to determine R I 0 , the basic reproduction number of system ( 8)-( 14) [30].The first step in using the next generation matrix approach in determining the reproduction number is identifying terms representing new infections in the model system and separating them from transfer terms [30].For this, let F i be the rate of appearance of new infections in compartment i, V − i be the rate of transfer of cells or parasites out of compartment i by all other means, and V + i be the rate of transfer of cells or parasites into compartment i by all other means.Set , and E 0 = X m f = (K, 0, 0, 0, 0, K i , 0) be the MFSS.Define The system is constructed such that x i = F i − V i , i = 1, 2, 3, 4, 5, 6, 7.Here X is the vector of our state variables.The matrix G := FV −1 is called the next generation matrix.The within-host basic reproduction number R I 0 is then largest or dominant eigenvalue of G.That is, R I 0 = ρ(G) = max λ {|λ | : Gx = λ x}.We recall here that it is not necessary to use the whole system to determine R I 0 , rather, only the subsystem that corresponds to the infectious compartments can be utilized [30].From our model ( 8)-( 14), the subsystem representing the infectious compartments are determined by the variables R p , M, G e , G l .Considering only the equations for these four variables R p , M, G e , G l , and identifying terms representing new infections in the system, separating them from the transfer terms yield Their corresponding Jacobian matrices, F and V , respectively, evaluated at the MFSS are Clearly, the matrix V is invertible.Through some standard computations we obtain where The expression for R I 0 in equation ( 16) is the basic reproduction number for the model system ( 8) - (14).We note that the same expression can be obtained by investigating the local stability of the MFSS.
In what follows, we will compute the normalized local sensitivity indices of R I 0 relative to the parameters associated to the within human-host malaria infection for model ( 8)-( 14).The set of input parameters for R I 0 is We will then investigate the impact on the normalized forward sensitivity (or elasticity index) of R I 0 due to the variations in the given parameters.When a system output has many parameters, the results of the output will not change equally to changes in parameters because some of them are highly sensitive, others have low sensitivity and some have relative zero sensitivity.The output can be optimized by examining which parameters are sensitive and which are not.
The normalized local sensitivity index of the output R I 0 , to a parameter p ∈ p, denoted here by Ψ .
Using this definition, we obtain the following indices for R I 0 relative to each of the parameters in (17): Remark 1. From the expressions of the local sensitivity indices above, and using the fact that 0 < σ < 1 we see that Thus, based on the sensitivity indices, a 1% increase in r, respectively β , while holding the other parameters fixed, will respectively produce a corresponding 1% increase in R I 0 .This is evident given the linear relationship between R I 0 and each of these parameters.A 1% increase in ω, respectively γ p , with all other parameters held fixed, correspondingly produces a less than 1% increase in R I 0 , meanwhile a 1% increase in each of µ p , μh , β 2 , µ p , ρ p , ρ m produces a corresponding less than 1% decrease in R I 0 .Hence, these are desirable target parameters for control, when it is feasible.
We are now in a position to use the computed normalized local sensitivity coefficients to investigates relative changes in R I 0 based on specified parameter values; see the example that follows.Example: For the following values, which were baseline values in [34], we obtained specified values of the normalized local sensitivity coefficients, summarized in Table 4 and illustrated in Figure 1.
The interpretation of the elasticity indices in the third column of Table 4 are as follows: For the specified parameter values given in (24), R I 0 is most sensitive to β 1 , the contact rate between HRBCs and free-floating merozoites, and r, the number of merozoites released per bursting IRBCs.They both have index value of +1.0, as expected from the computation of the index function.Variation or uncertainty in these parameters will bring more variations or uncertainties in R I 0 and consequently, the outputs of the model.In particular, for all variables, if 10% more (respectively less) merozoites are released (that is r changes by 10%) or come in contact with HRBCs (that is β 1 changing by 10%), R I 0 increases (respectively decreases) by 10%, each.The next set of parameters to which R I 0 is most sensitive to are: ω = Λ − µ h , the per capita linear growth rate; μh , the death rate of HRBCs due to density dependent contact inhibition, and µ m , the per capita natural death rate of merozoites.Their respective normalized elasticity indices are ≈ +0.9505, − 0.9505 and −0.9498.Among the three parameters, the one most feasible to control is µ m , with index value of ≈ −0.9498.Thus, with all other parameters held fixed, a 10% increase (respectively decrease) in the natural death rate of free-floating merozoites, produces a corresponding decrease (respectively, increase) of approximately 9.5% in R I 0 .The index values for the other parameters, given in order of decreasing magnitude are: γ p , rate at which IRBCS change paths, either bursting to release merozoites or continuing towards the gametocytogenesis path; K i , the carrying capacity of innate immune cells; ρ p , the elimination rate of IRBCs by innate immune cells; β 2 , the absorption rate of free merozoites by IRBCs; µ p , the elimination rate of parasitized red blood cells; σ , the proportion of merozoites that continue the path towards gametocytogenesis and ρ m , the elimination rate of merozoites by innate immune cells.They each have magnitude less than 0.1 as shown in Table 4 and can be described similarly.What is clear is that if we focus on parameters that can be controlled in order to reduce R I 0 and thereby reducing malaria parasitemia within a human-host, we need to reduce the infection rate β 1 , the number of merozoites released per bursting IRBC or increase the per capita natural death rate of free merozoites, µ m .

Graphical results showing the local impact of the parameters
on R I 0 First, we present a plot of the sensitivity indices of the parameters that are associated to R I 0 .See Figure 1.The model parameters influence the size of the reproduction number R I 0 in different ways.Moreover, the sensitivity indices shown in Table 4 for the parameters whose indices differ from 1, will change when some of the other associated parameters in the computed index change.Thus, it is important to investigate how the change in R I 0 is manifested with an increase or decrease in a specified parameter, or a combination of parameters.Beginning with µ m , while keeping other parameters fixed, we plot R I 0 and also the size of the normalized local sensitivity coefficients Ψ R I 0 p against µ m .The results are illustrated in Figures 2 and 3 for the case where the parameters ω, μh , γ p , r, σ , K i , ρ p , ρ m are as stated in the parameter example equation (24), while β 1 and β 2 are each increased 10 folds so that β 1 = 5 × 10 −6 = β 2 and µ p is increased from 0.0091 to 0.07.From the plots, as the free merozoites per capita death rate, µ m , is increased from 24 to its maximum value of 72, keeping 96|.This should be expected because a higher µ m corresponds to a shorter window for parasites to infect HRBCs.This illustrates that we can achieve parasite clearance by controlling µ m .This can be done via drugs that speed up the death rate of free merozoites within the human blood stream.increases, however, its value decreases with increase in µ m , which confirms the previous results.So one of the ways to eliminate or reduce the parasite density within the human blood is to shorten the life expectancy of merozoites so that they can die before they have the opportunity to infect HRBCs.Another way is by inhibiting successful contacts between merozoites and HRBCs.Similar plots can be obtained for the other parameters to see their effect on reducing the reproduction number R I 0 .From the normalized local sensitivity indices computed in equations ( 18)- (22), it is evident that although these indices lie in [−1, 1], their magnitude, whether closer to one or closer to zero will depend on the other parameters present in the expression.
As an example, Ψ R I 0 γ p = ρ p K i +µ p ρ p K i +γ p +µ p as shown in equation (20), and clearly depends on ρ p K i + µ p .As ρ p K i + µ p → 0 (which can occur for small ρ p and µ p values), Ψ R I 0 γ p → 0. If, on the other hand, ρ p K i + µ p → L → ∞, (which can occur for large K i values and reasonable ρ p and µ p values, Ψ R I 0 γ p → 1.We note here that the statement → ∞ just represents, from a mathematical perspective, a large value within feasible biological ranges.We illustrate this point by plotting 3D surfaces of some of the µ m (ρ p , ρ m ) (4a) and R I 0 (ρ p , ρ m ) (4b) as functions of ρ p and ρ m , with ρ p ranging from 0.78 × 10 −6 to 0.98 × 10 −3 and ρ m ranging from 0.9 × 10 −5 to 0.9 × 10 −1 .All other parameters are as given in equation (24).Fig. 5: Contour plot of R I 0 (µ m , β 1 ) as functions of ρ p and ρ , with µ m ranging from 48 to 72 and β 1 from 5 × 10 −9 to 3 × 10 −6 .Except for β 2 = 5 × 10 −6 , all other parameters are as given in equation (24) elasticity indices as functions of two input parameters that co-impact the index; see Figures 6-9.
In Figures 7-9 we have shown that in general, when two parameters are combined, their combined effect on R I 0 is generally stronger and will have a stronger effect on impacting parasitemia.

Discussion and conclusion
We set out to investigate the local sensitivity analysis of the threshold disease parameter, the basic reproduction number R I 0 , for a within-human host model of the immuno-pathogenesis of the malaria parasites, originally developed in [34].By computing the normalized local sensitivity coefficients for R l 0 , we investigated the impacts of certain key control parameters and the relative impacts of pairs of parameters to this threshold value when other parameters are held fixed.We also gave a brief review of two local sensitivity methods: the indirect finite-difference method and the direct differential method.The direct differential method is a fundamental local sensitivity analysis method which involves calculating partial derivatives of the model output with respect to the input parameters.and defined in equation (20), with respect to β 2 and ρ m (graph (a)), and also with respect to ω and ρ m (graph (b)).In Sub Figure 7a, β 2 ranges from 5 × 10 −7 to 5 × 10 −6 and ρ m ranges from 10 −8 to 0.1, with all other parameters as in (24).In Sub Figure 7b, ω ranges from 0.036 to 0.5, with β 2 = 5 × 10 −6 and all other parameters fixed as in (24).In general, sensitivity analysis allows for, when possible, the exact computation of the sensitivities of a systems response to variations in the systems parameters, around their nominal values, which then allows for a discussion on the relative influence of the input parameters to the model results.With this, a discussion on parameters that are less or more important can be presented which could provide information on whether a parameter may be eliminated from the model system [6,12].However, it is worth noting that even through a parameter may have great influence on a system's output, the question of whether in reality that parameter may be feasible, easy, efficient and cost effective to control is worth investigating.
For the within-human-host immuno-pathogenenic malaria model presented, we applied the direct differential method to investigate the local sensitivity of R I 0 to the model parameters.The parameters affect the densities of HRBCs, IRBCs, free merozoites, the early/immature stage (stages I -IV) and late/matured stage (stage V) gametocytes, as well as the functioning of the human innate immune and adaptive immune system effectors.The basic reproduction number, R I 0 , was computed using the next generation matrix (which matched the threshold parameter that determines the existence of the parasitized steady state solutions), from which the normalized local sensitivity coefficients in relation to each parameter were computed.For a representative feasible parameter set, the numerical values of the normalized sensitivity coefficients were computed as summarized in Table 4.
Based on the computed indices, we can explore the impact of varying a parameter value on the ability of the parasite to persist in a human, that is when the reproduction number R I 0 surpasses 1.This is important as it can inform on possible parasite control parameters since normalized local sensitivity coefficients quantify the percent change that results in an output variable based on a one percent change in an input parameter, keeping other parameters fixed [29,38].Note that these normalized local sensitivity coefficients are dimensionless and highlight relative changes.For our model, the normalized indices are shown on Table 4. Negative normalized local sensitivity coefficients indicate that R I 0 is indirectly proportional to the corresponding parameter and so an increase (respectively decrease) in the parameter would produce a corresponding decrease (respectively increase) in R I 0 , of size relative to the size of the index value.On the other hand, parameters with positive index are directly proportional to R I 0 , and changes are in the same direction so that an increase (respectively decrease) in the parameter would produce a corresponding increase (respectively decrease) in R I 0 .In our model, parameters that are biologically feasible and reasonable to control are the free-floating merozoite death rate µ m , which has a large negative elasticity index.Thus with other parameters fixed, speeding up the deaths of free merozoites will have a significant impact on the reduction of R I 0 .More specifically, a 10% increase in µ m produced a corresponding 9.5% reduction in R I 0 when other parameters are as shown in (24).From Figure 8b, the size of this impact is muted and can be enhanced as well when the sizes of other parameters are considered.Another parameter of interest is γ p , which determines the bursting time frame for parasitized red blood cells and hence successful passage of parasites towards the sexual stages, the gametocytes, for those IRBCs that do not burst.From Table 4 the index value for γ p is 0.0746 indicating that a 10% reduction in this parameter produces a corresponding 0.75% reduction in R I 0 when other parameters are as shown in (24).Both of these parameters are biologically feasible control parameters.
To highlight the previous discussion, the Plasmodium falciparum parasites exhibit different life cycle stages.Stopping or slowing one of these stages will stop or slow down the production of merozoites within a human or the development in a human as well as its transmission to mosquitoes.This is mostly achieved via antimalarial drugs.However, in recent years, significant increase in drug resistance has propelled researchers to expand their research questions and endeavor and consider other strategies that can either speed up merozoites killing, slow gametocyte development and persistence in the bloodstream, as well as find ways to block transmission to mosquitoes.These control measures can occur through the use of effective non-resistant anti-malarial drugs or vaccines (various vaccine constructs are undergoing evaluation in clinical trials while others are in advanced preclinical development phase [33]) designed to disrupt parasite reproduction and further development in the mosquito's midgut, and can thus serve to break the life cycle of the parasite [13].A vaccine or effective drug that kills free merozoites quickly (an increase in µ m by shortening the life expectancy of merozoites) or seeks to inhibit bursting of IRBCs (a decrease in γ p ) or seeks to inhibit the ability of merozoites to commit to switching to the sexual stages or one that inhibits the successful maturation of immature gametocytes within a human or a successful transmission of gametocytes from human to mosquito will be desirable and can play a great role in reducing and possibly eliminating parasitemia in a human.This would produce a corresponding positive effect in reducing gametocyte transmission to mosquitoes.Hence a study as the one done here, is crucial in identifying important and feasible control parameters.
Our model did take into consideration the behaviour of the innate and adaptive immune behaviors.Feasible and realistic parameters that can activate the immune cells and enhance their effect in killing parasites are important.From Table 4, the elasticity index value for ρ m is small indicating that its individual impact on R I 0 may be small.However, the influence of µ m on R I 0 can be dampened or enhanced when both β 2 and ρ m (see Figure 6b) parameters are varied.Hence, a combine control strategy that inhibits contact between HRBCs and merozoites as well as increases the efficiency of innate and adaptive immune responses in killing, blocking, inhibiting parasites growth would be desirable.A way this can potentially be achieved is through vaccines that can enhance the activities of the innate and adaptive immune system in their actions against parasites.
From our discussions above our normalized local sensitivity analysis for R I 0 does give insight to the significance of some parameters in reducing malaria parasitemia through their effects on R I 0 .It also highlights which parameters are highly correlated with the output R I 0 .We comment that although there are generally two broad techniques of sensitivity analysis, local sensitivity analysis and global sensitivity analysis, we elected to use a local sensitivity analysis in our work.As noted in [6], an objective of local sensitivity analysis is in the analysis of local behaviours of a system's response around a chosen point or trajectory in the combined phase space of parameters and state variables.Thus it focusses on the impact of small perturbations on the model outputs [6].On the other hand, global sensitivity analysis helps to understand how the model outputs are affected by large variations of the model input parameters when all the input factors are varied simultaneously.The sensitivity is evaluated over the entire range of each input factor [38].We believe that a local sensitivity analysis for the within-human host parasite dynamics is a reasonable choice because one can see the impacts of small variations.Moreover, given that the parasite forms are very different and behave differently in their residing systems, a local sensitivity analysis provides reasonable insights as to changes in parameters to the model outputs within each local system.See [19,29] for some examples on the use of local sensitivity analysis (LSA) in biological and physical systems.
Fig. 7: A 3D plot of the local sensitivity index function of R I 0 to β 2 , denoted by Ψ R I 0 β 2 Fig. 8: A 3D plot of the local sensitivity index function of R I 0 to γ p , Ψ R I 0 γ p defined in equation (20).Sub Figure 8a is a plot with respect to γ p ranging in [0, 1] and ρ p ranging in [10 −8 , 0.1] while Sub Figure 8b is a plot with respect to γ p ∈ [0, 1] and µ p ∈ [0.036, 0.2].All parameters are held fixed as in as (24) with β 2 = 5 × 10 −6 .

Table 1 :
Summary comparison of Local sensitivity vs Global sensitivity

Table 2 :
Description of state variables and their quasi-dimension: C = the density of red blood cells per unit volume of blood =red blood cells /µl; M = density of free-floating merozoites per unit volume of blood = merozoites/µl; G = density of gametocytes per unit volume of blood = gametocytes/µl; I = density of immune cells per unit volume = immune cells /µl.

Table 3 :
Description of Parameters and their quasi -dimensional units

Table 4 :
(24)l sensitivity indices of the reproduction number for the general model with immunity uysing the parameter values given in(24).