The sickle cell trait originates as a natural mutation of the hemoglobin gene. Such mutation results in replacement of glutamic acid at position 6 of the beta chain (β6) of hemoglobin by valine. The mutation results in the aggregation, in the form of a polymer, of the sickle cell hemoglobin (HbS) when it is in the deoxygenated state [1]. The polymerization process takes place in two stages, which are separated by a time delay [2]: 1) homogeneous polymerization where monomers join to form a polymer; 2) heterogeneous polymerization where monomers join an existing polymer. The polymers formed in the first step can melt yielding monomers if they return to the lungs quickly enough to be re-oxygenated. The polymers that do not return to the lungs in a timely manner will likely go through the second stage. The mechanisms for the homo-and heterogeneous polymerization are referred to as single and double nucleation [3].
Ferrone's work [3] focused on the growth of sickle hemoglobin (HbS) polymers or fibers. He explained that the growth of the HbS fibers follows the double nucleation mechanism, which he described as homogeneous and heterogeneous nucleation. The formation of a polymer is initiated by homogeneous nucleation in the solution phase. Local temperature and concentration of the homogenous nuclei can then increase the formation of additional polymers on the surface of an existing polymer (heterogeneous nucleation). Local temperature is very relevant to our work because we are looking at the change in the rate constant where temperature is the only factor that causes the change in parameters. The surface area of heterogeneous nuclei is constantly increasing with time [3]. From the understanding of these two processes, Ferrone developed a mathematical model for the melting of the polymers (de-polymerization) on the assumption that the polymer melting is the reversal of polymerization. The model can be described by two rate equations, one for the formation of polymers and the second for the incorporation of monomers into polymers.
The kinetics of the sickle cell hemoglobin polymer melting has been studied experimentally by using the stopped flow method where the melting was monitored using light scattering [2]. The results showed that polymers melt more quickly in cells containing oxygen or carbon monoxide. Therefore, two sets of experiments were conducted to study polymer melting. The first experiment involved monitoring the rate of HbS polymer melting in a deoxygenated phosphate buffer at pH 7.1, 25 • C, and in the second experiment the same buffer was saturated with carbon monoxide (CO) [2]. The author concluded that polymers melt more efficiently in presence of CO. Additionally, the authors noted that for the model to fit the experimental data, CO was assumed to bind to HbS in the polymeric and in the monomeric form [2]. In Aroutiounian's study [2], the model was rewritten as a system of two equations by including the fact that polymer only melts from the ends. Thus, the study improved Ferrone's model, which is based on a single differential equation. Further investigation of the effects of CO binding to the depolymerization of solution phase of the polymer used a nonlinear model that extended the original two-species model [2] into a model that included the four species: deoxy HbS monomers, CO-bound HbS monomers, deoxy Hbs polymer and CO-bound HbS polymer [4]. The model assumed that CO binds to both the monomeric and polymeric forms of HbS and described the dynamic interaction in each phase of the melting and their CO binding processes. Results from the analyses predicted that melting of de-oxy polymers occurs rapidly when the solution was saturated with CO. Additionally, they observed from the model equations that not all the polymers melt in the CO-binding equilibrium stage, indicating that the melting of the CO-bound polymer takes place as an equilibrium process.
Normally cell behaviors are determined by the interaction of the components in the biological system instead of the characteristics of the individuals. The impact on the system output from various parameters is usually different. In biological experiments, the impacts from different parameters are often difficult to determine as it requires the repeated experiments and the subtle measurements which may be practically impossible [5]. This, however, can be easily achieved in mathematical modeling with numerical simulations by parameter sensitivity analysis. Sensitivity analysis is the study of how the uncertainty in the model output can be attributed to different sources of uncertainty in its inputs [5,6]. This sensitivity analysis provides guidance toward the parameters that must be taken into consideration during the experimental design [7].
There have been studies concerning parameter estimation, which have led to approximation of parameters value in the nonlinear dynamic system [2,4]. However there remains the possibility to improve the measurement of the parameter values. In this paper the importance of parameters in the extended model of melting (de-polymerization) of HbS polymers in the presence of CO is assessed.
There are many approaches to study the sensitivity of a model output with respect to the input parameter [5]. A comprehensive evaluation of various sensitivity analysis methods (including the FAST and Sobol methods) is provided in a case study on a hydrological model [8]. This study particularly focuses on a traditional sensitivity analysis (TSA) and a multi-parameter sensitivity analysis (MPSA) on the extended version of the model with four species [4]. All four concentrations as model outputs and four binding/melting rates as input parameters are considered in this study. The sensitivity analysis results from TSA and MPSA are essentially consistent.
The remainder of this paper consists of the following sections. Section 2 presents the derivation of the mathematical model equations starting from a description of the double nucleation process to the extended CO-mediated sickle cell polymer melting (depolymerization). Then, Section 3 explains the detailed methodology of the parameter sensitivity analysis including the TSA and MPSA. Section 4 reports and discusses the obtained results from our analysis. Section 5 summarizes the conclusion of the study and provides the recommendation for future work. Lastly, Appendix includes the detailed sensitivity equations for TSA.
Our current version of the model is an extension of the basic model proposed in [2,3]. The model developed for the HbS polymer melting is based on the observation that the HbS fiber melting is the reversal of growth. According to [3], the growth of fiber occurs through a double nucleation process. Figure 1 depicts the process.
As it happens in the solution, homogeneous nucleation is the simple process where the monomers attach to each other to form a polymer. In heterogeneous nucleation, the monomers add onto an existing polymer to form a more complex polymer. As Fig. 1 depicts, the length of the arrows going back and forth changes in the different stages of the figure, which indicates the change in the relative rates of fibers melting or growth. In the early stages of the homogeneous nucleation phase the rate constant of melting (K backward ) is higher than the rate constant of growth (K forward ). The equilibrium stage (K eq ), where the rate constants of the growth and melting are the same (K forward = K backward ), happens when the polymer reaches a certain size. Beyond this size, the rate of the polymer growth becomes higher than the rate of melting. The driving force for this switch in kinetics is the increase in the stability of the polymer as it exceeds the above-mentioned size. This increase in stability is due to the increasing amount of energy released as more and more monomers join the polymers which makes the polymer formation to have higher negative ∆G (meaning polymer formation becomes thermodynamically more favorable). The same cascade of events takes place in the heterogeneous nucleation phase. Again, when the size of the polymer exceeds a certain limit the formation -or growth-of the polymer becomes more favorable than the melting of the polymer. The final size of the polymer is limited by its solubility. This means once the size of the polymer -the fiber-in the heterogeneous nucleation exceeds a certain limit it becomes insoluble, or a better explanation is that the polymer stops growing when there are no more monomers to bind to it. When the size of the polymer reaches this limit it would deplete the solution from any monomers -now the concentration of the monomer to begin with is limited by its solubility at the conditions at hand (temperature, phosphate buffer concentration and the fact that the buffer is saturated with carbon monoxide).
The double nucleation model can be described by two rate equations, one for the homogeneous and the other for the heterogeneous. Since in either homogeneous or heterogeneous nucleation there is an addiction of a monomer to the nucleus, they both can be represented as (this applies to the CO bound and deoxygenated polymer) the equilibrium P i + (monomers) ⇐⇒ P j .
Here P i and P j are two forms of the polymer where the j form has one monomer more than the i form, or P i is the homogeneous polymer and P j is the heterogeneous polymer. The formation of the polymer is thus expressed by:
and the rate of disappearance of the monomers from the solution phase into polymer is given by:
where C p,i (t) and C p,j (t) are the time-dependent concentration of polymers i and j, respectively, and C m (t) is the monomer concentration, k + is the concentration-independent rate constant for addition of monomers to nuclei or polymer i , and k -is the concentration-independent rate of the dissociation of a monomer from polymer j.
Briefl [9] observed that the growth is an elongation and the melting is a shortening of the HbS fiber at the ends. Thus melting of the HbS polymers can only occur at the end of the polymers. On a short time scale, the concentration of polymers P i and P j remains constant at a certain time. Therefore, Aroutiounian [2] sets the concentration of polymers P i and P j to be the same as C p (t); he replaced k + with k- Cs , where C s is the solubility concentration of the de-oxygnated HbS polymer in the buffer at T , and k -with k d indicating the rate constant of the melting of the de-oxygenated HbS polymer. So Eq.2 is rewritten as the following in term of the melting rate:
Also, the total HbS molar concentration (C tot ) is the sum of the molar concentration of the hemoglobin molecules in the polymer phase, C p , and that in the monomer phase, C m . C tot (t) = C m (t) + C p (t). After differentiating and substituting into Eq.3 with the fact that dCtot(t) dt = 0, we have:
C. The Extended Model with Four Equations With the CO-mediating the melting of the polymer incorporated [4], the polymerized and monomerized populations are then divided into two sub-populations: C CO p (t) and C CO m (t) which are CO-bound polymer and monomer HbS, respectively; and C d p (t) and C d m (t) which are de-oxy polymer and monomer HbS. Since CO binds to the monomer and polymer tightly, the model then assumes that CO-binding results in a decrease from de-oxy HbS C d m (t) in the solution phase, which then becomes a gain for the CO-bound solution phase C CO m (t), with a CO binding rate constant k m . Taking all these into consideration produces Eq. 5. Likewise with the polymer phase molecules, the assumed CO-binding outcome is a loss from the C d p (t) and a gain to C CO p (t) by the CO binding rate constant k p . Equation 6 corresponds to the production of COmonomer. It is generated while making the same assumptions for the melting of the CO-polymer as the ones made when writing the differential equation of the melting of de-oxypolymer (Eq. 3), with the melting rate constant being K CO , and the formation of the COmonomer from deoxy monomer in the presence of CO with the rate constant being K m . Equations 7 and 8 are produced the same way as Eq. 4 with taking into consideration the binding of CO to deoxy polymer. Our model is then based on the following:
These models allow CO-binding to polymers and melting occur at the endpoints as well as at the surfaces of polymer fibers [4]. The diagram in Fig. 2 describes the reaction paths in the model equations. With CObinding, the de-oxygenated polymers/monomers of HbS produce CO polymers/monomers HbS. The melting/depolymerization of the CO-bound or de-oxygenated polymers yields CO-bound de-oxygenated monomers, respectively in an equilibrium process.
To mathematically analyze the system, we simplify the notations by replacing ( [10]. We then have the following model equations:
III. METHODOLOGY
The traditional sensitivity functions (TSFs) are used in the traditional/classical sensitivity analysis. To investigate and quantify the sensitivity of a model output resulting from the variations in the parameters, we consider the partial derivatives of the output variable with respect to the parameters that it depends on [11,12]. The detailed TSFs for the sickle cell model in this study are derived and presented as follows.
Many times, the model is described by a system of ordinary differential equations (ODEs). The system can be written in the vector form:
where X ∈ R n denotes the state variable vector and P ∈ R r denotes the vector of parameters and the initial condition is X(t o ) = X o . Depending on the research goal, the output variable vector Y for the sensitivity analysis can be a subset of X or the full state X. By sensitivity we mean how Y changes with respect to P , that is ∂Y (t) ∂P . Since all the states are related and coupled, solving ∂Y (t) ∂P requires solving the full states
∂P . In the following, we work on the full state vector as the output variable vector. Note that as X is a function of time, the sensitivity is also a function of time. By differentiating both sides of Eq. 13 with respect to P , we have a system of differential equations for the sensitivities:
Now if we reverse the order of differentiation and couple the equation with Eq. 13 then we get a n + nr dimensional system of ordinary differential equations for both the model variables and the sensitivities dX(t) dt = F (X, P )
Here, we assume that ∂X(0) ∂P = 0, because the initial conditions for the model would be considered independent of the parameters.
We apply the above Eqs. 15 and 16 onto our sickle cell model Eqs. 5-8. The new system is obtained with 4 original model equations and 16 equations for the TSFs. The detailed equations for the TSFs are presented in Appendix.
First developed in the field of hydrology [13][14][15], the MPSA method is a technique that quantifies the relative importance of the parameters related to the output variables in the model [16]. The MPSA considers the combined effects of multiple parameters on the output of the model. We first use the Monte-Carlo method to randomly select values from the distributions of the considered parameters. Similar to the studies reported in [17,18], we use the uniform probability distribution since the natural distribution of the parameters in the sickle cell polymer decomposition system is unknown. With the randomly selected parameters values, the differential equations model is then simulated repeatedly. Next, an objective value is calculated by using the objective function (defined below) to classify the output of the model simulation as either acceptable or unacceptable. Lastly, a statistical evaluation is carried out for the acceptable and unacceptable cases, giving quantified values for the sensitivities of the parameters.
The detailed procedure of the MPSA can be found in Cho et al. [16]. In this study, instead of setting the range between one fifth of a nominal value and five times the nominal value [16], for the sickle cell model, we set the range between one half of the nominal value and two times the nominal value since roughly speaking the rates may be doubled or halved when the temperature changes within ten degrees. In comparing the two cumulative distribution functions (CDF) of the parameter values associated with the acceptable and the unacceptable results, Cho et al. [16] simply consider the "cumulative frequencies" for each parameter via corresponding correlation coefficients. In this study, a more thorough comparison technique is adopted, namely, the Kolmogorov-Smirnov (KS) test [16]. The KS distance is calculated, with the large distance value indicating the large sensitivity of the parameter, since the large KS distance implies that the two CDFs are different to each other.
The objective value is defined as the sum of squared differences between the output values from the sampling parameters and the output values from the nominal parameters:
where f obj (k) is the objective function that describes how much the system output with the sampling parameters changes from the data with the nominal parameters; f nominal (i) denotes an output value from the nominal parameters at the ith time; f sampling (i, k) denotes the output value from the sampling parameter k at the ith time; and q is the number of time point.
A typical solution behavior of the model over time and the general dynamics of the CO-mediated depolymerization process can be seen in Fig. 3. The de-oxy monomer x(t) first increases to its peak, then decreases gradually to 0, while the de-oxy polymer z(t) decreases the whole time to 0. The CO-bound monomer y(t) first increases to its peak which is higher than the solubility, then decreases to its solubility, while the CObound polymer u(t) increases the whole time to its solubility. For the solutions in Fig. 3, the parameters are set as k 1 = 0.028, k 2 = 0.07, k 3 = 0.1, k 4 = 0.01. The solubilities are C 1 = 0.4 and C 2 = 0.8. The initial conditions of the system are chosen as x(0) = 0.0036mM, y(0) = 0, z(0) = 1.175mM, u(0) = 0. A detailed analysis with rigorous mathematical proofs on the dynamics of the extended model Eqs. 9-12 can be found in Daniels-Jones et. al [10].
The above parameter values are set as the nominal values of the parameters. The above solubilities and initial conditions for Fig. 3 The initial values and parameters ranges are chosen based on previous experimental studies on the COmediated sickle cell de-polymerization [2,3]. Among all species, the initial concentration of de-oxy polymers z(0) is the highest (1.175mM), because medically this is the species that causes sickle cell illness. The initial concentration of de-oxy monomer x(0) being low value (0.0036mM) is justified due to the high value of de-oxy polymers. The initial concentrations of the other two species (y(0) and u(0)) are assumed to be 0 because they are not exposed to the CO yet. Based on the model equations, the results remain similar for some range of initial and parameters values around the chosen ones. For very different initial and parameters values, the results of sensitivity may be different, which is not covered in this study. For the numerical simulations, we use the Runge-Kutta method for the time-marching integration of the ODEs, with time step dt = 0.01 and the time interval [0, 400], which provides an adequate accuracy and sufficient transient time.
People often vary the parameter values and conduct experiments for responses to observe the effects of the parameters over the output variables. A rough estimate and some findings may be obtained in this way. As an example, with a few different k 1 values (0.0112, 0.0224, 0.0336, 0.0448, 0.0560), the corresponding solutions of all four concentrations are presented in Fig. 4. From Fig. 4a, the increase of k 1 leads to the increase of the peak for x; the peak time remains the same; after the peak, the decreasing speed is large with a large k 1 -value. From Fig. 4b, the larger the k 1 -value is, the y-value increases faster to its higher peak (over the solubility 0.8), then decreases slower to 0.8 the solubility. From Fig. 4c, the larger the k 1 -value is, the faster/quicker the z-value decreases to 0. Lastly, the u-value first increases slower with larger k 1 -value. At some point of time, this is reversed: the u-value increases faster with larger k 1 -value, as shown in Fig. 4d. It takes shorter time for the u-value to reach its solubility 0.4 with a larger k 1 -value. The rough effect of the k 1 -value on the individual variable can be observed in Fig. 4. However, some detailed sensitivity analysis, such as on which variable does the parameter k 1 have the most effect, cannot be obtained from Fig. 4 directly.
The effects of all four parameters over the variable CO-bound polymers concentration u are displayed in Fig. 5. All four parameters have some effects on variable u. From Fig. 5a and 5c for parameters k 1 and k 3 , it starts with the lower the parameter values, the faster the u-value increases; after some time, this is reversed: the lower the parameter values, the slower the u-value increases. The solutions with various parameter values for k 2 and k 4 in Fig. 5b and 5d show the consistent trend: the lower the parameter values, the slower the u-value increases. The rough effect of the individual parameter over the u variable can be observed in Fig. 5. However, some detailed sensitivity analysis, such as which parameter has the most effect on the variable, cannot be obtained from Fig. 5 directly.
For a systematic sensitivity analysis and to compare the impacts of the parameters, we utilize the TSA and the MPSA to examine the sensitivities of the output variables (x, y, z, u) that are presented in our model equations with respect to the changes in the parameters' values (k 1 , k 2 , k 3 , k 4 ). In the following 3 subsections, we quantify the sensitivity with estimates, and we rank the effects of parameters accordingly.
A. Results from TSA 1) TSFs with parameters at nominal values: With the nominal values for the parameters, the ODEs for the TSF functions in Appendix can be solved numerically. Figure 6 shows the TSF functions ∂x ∂k1 , ∂y ∂k1 , ∂z ∂k1 , ∂u ∂k1 of the variables with respect to parameter k 1 . The behavior of the TSFs in Fig. 6 is consistent with the phenomenon in Fig. 4. The larger the TSF value is, the larger the sensitivity (change/variation) of the variable is. The positive TSF value means that the increase of the parameter leads to the increase of the variable value, while the negative TSF value means the opposite. From Fig. 6a for the sensitivity of variable x over parameter k 1 , the sensitivity increases to its peak at t ≈ 20, then decreases to 0 at t ≈ 75, continues to decrease and then increases but remains negative, lastly at t ≈ 200 it settles at 0. The time stamps (20, 75, 200) match well with those in Fig. 4a. Similar consistency can be found in Fig. 6b vs Fig. 4b, in Fig. 6c vs Fig. 4c, and in Fig. 6d vs Fig. 4d. It should be noted that all the sensitivities converge asymptotically to 0 over some time. Note that instead of the discrete values for the parameter in the trials in Fig. 4, the TSF describes the instantaneous change/impact of the variable with respect to the change (increase or decrease) in the parameter. Therefore, as a quantified sensitivity, the TSF in Fig. 6 is more rigorous and more reliable than the rough estimate observed from Fig. 4. The sensitivities presented in Fig. 6 are functions of time, i.e., the sensitivity values are different at different times. To better compute the sensitivity of the variable to the parameter, we sum up the sensitivity values over the time. For this analysis, we use the L 2 norm of the TSF to measure the sensitivity. Note this norm is defined as [11]:
Since we are approximating the solution to our sensitivity ∂X ∂P , we also approximate this norm. The L 2 norms of the TSFs are shown in Table I.
With the quantified sensitivity values in Table I, we can compare and rank the sensitivities of variables to the parameter by reading the table in columns. The column of k 1 for the sensitivities of variables to the parameter k 1 shows that variable z is the most sensitive, followed by y. The sensitivities of u and x variables are similar. The sensitivity values in column k 1 quantify the phenomenon in Fig. 4, sum up the values in Fig. 6, and are consistent with the figures. Note that the scales in the vertical axes in Figs. 4 and 6 vary. All the sensitivity values in column k 2 are smaller than those in column k 1 , which implies that k 2 may not affect the variables as much as k 1 . The effect of k 3 is simple with 0 sensitivity for variables x and z, and some small sensitivity for the other two variables. The values in k 4 column are large, indicating the large effect on all variables. In a summary, from the columns in Table I, the order of variables' sensitivities (from large to small) is (z, y, u, x) for k 1 ; (x, y, u, z) for k 2 ; (u -y, x -z) for k 3 ; and (u, z, y, x) for k 4 , where the dash "-" means the same sensitivity. ∂k 1 ; b) the TSF of variable y(t) with respect to k1, ∂y(t) ∂k 1 ; c) the TSF of variable z(t) with respect to k1, ∂z(t) ∂k 1 ; d) the TSF of variable u(t) with respect to k1, ∂u(t) ∂k 1 .
Reading in rows of Table I, we can compare and rank the sensitivities of the variable to parameters. In row x, parameter k 4 is with the largest sensitivity, followed by k 1 and then k 2 . k 3 has no effect on variable x. From the other rows for the other variables, even though the sensitivity values are different, but the order of parameters is consistent: k 4 , k 1 , k 2 , k 3 . Note that the values in row u in Table I quantify the variation in Fig. 5 for variable u with different varying parameters. The measure provided in Table I is more reliable and clearer than the rough variation in Fig. 5. The order/rank of the parameters (k 4 , k 1 , k 2 , k 3 ) may not be obtained directly from Fig. 5.
2) TSFs when all parameter values vary: The above sensitivities are measured for the nominal values of the parameters. With different parameters' values, the changes on the variables are different. Figure 7 shows the variations of the variables for various parameter values in the ranges around the nominal values. The reference variable values are the solutions when the parameters are at the nominal values. For Fig. 7, the parameters' values are chosen randomly from the ranges of the parameters. Therefore, the effects of the combinations of the parameters are presented in Fig. 7.
It can be seen from Fig. 7 that the solution values vary above and below the reference values with different values for the parameters. However, it is not easy to analyze the sensitivity behavior directly from
Variables k1 k2 k3 k4 x 41.5601 25.9384 0 44.3415 y 77.1920 20.7210 10.3865 100.7236 z 113.7541 1.9437 0 128.8967 u 47.7038 20.0600 10.3865 203.1208 Table I: TSA results: the sensitivities when parameters are nominal values. Variables k1 k2 k3 k4 x 0.6195 0.9640 0 0.0772 y 1.0130 0.8312 0.4162 0.7778 z 1.5828 0.0025 0 0.7172 u 1.0570 0.9943 0.4162 1.9409 Table II: TSA results: the sensitivities when all parameters vary.
Fig. 7. To quantify the measure of the sensitivities, we utilize the TSFs. As an example, Fig. 8 shows the TSFs of all variables with respect to parameter k 1 when all parameters are varied. The behavior of the graphs follows a similar pattern due to the increase or the decrease of the parameters' values.
Next, to wrap up the results in Fig. 8, we seminormalize the TSFs. Normalization is a process that is used to eliminate redundancy, reduce the potential for anomalies during data processing and maintain the consistency and integrity of the data. The L 2 norm of the TSF with the sampling parameters referenced to the TSF with the nominal parameters is given
• ∆t where X = {x, y, z, u} , P = {k 1 , k 2 , k 3 , k 4 } , P s is the sampled parameter value, and P o is the nominal value. The semi-normalizer consists of multiplying the norm by the distance between the parameters. It is defined as follows:
The results, the semi-normalization norm for the TSFs of each variable and each parameter is reported in Table II. The values in Table II measure the difference between the magnitudes of TSFs. The larger the value is, the larger difference between the magnitudes of the TSFs is, thus the more sensitive the variable is to the parameter. Fig. 8: The TSFs for all variables with respect to parameter k1 when all parameters vary. (a) the TSF of variable x(t) with respect to k1, ∂x(t) ∂k 1 ; b) the TSF of variable y(t) with respect to k1, ∂y(t) ∂k 1 ; c) the TSF of variable z(t) with respect to k1, ∂z(t) ∂k 1 ; d) the TSF of variable u(t) with respect to k1, ∂u(t) ∂k 1 .
k 1 for variable y are presented in the Fig. 9, and their CDFs of acceptance and non-acceptance are displayed in Fig. 10. The PMFs and CDFs of other parameters are similar to Figs. 9 and 10. The distributions are tested, and the Kolmogorov-Smirnov distance is computed to rank the sensitivity of each parameter, with the results shown in Table III. From these results including the figures and the table, for the output variable y the ranking of all four parameters sensitivities is
From the values in Table III and the graphs in Figs. 11 and 12 for the PMFs and CDFs of all four parameters for variable y, the sensitivity of parameters k 3 and k 4 are very close with no significant difference.
Similarly to the analysis on the output variable y in the above, we compute the PMFs of acceptance and non-acceptance, the CDFs of acceptance and non-acceptance, and the Kolmogorov-Smirnov distances for each output variables x, z, and u. Since the rank of the sensitivity of each parameter remains the same in each chosen percentile (the 25 th , the 37.5 th , the 50 th , the 62.5 th , the 75 th ) and the mean of the objectives for one variable, we choose the mean values of the objectives for each variable to report in detail in Table IV. From Table IV, with the output variable y the ranking for all four parameters sensitivities is
In comparison of the results from the TSA and MPSA methods, Table V displays the ordering of the ranks. Here, we notice that both results are consistent except for the output variable y whose 3 nd and 4 rd sensitivity ranks are reversed in the TSA and MPSA methods, and the variable u whose 2 nd and 3 rd sensitivity ranks are also reversed in the TSA and MPSA methods. Nevertheless, the variable x is most sensitive to changes in the parameter k 2 , the variable y and z are most sensitive to changes in the parameter k 1 , and the variable u is most sensitive to changes in the parameter k 4 . k 3 does not affect any output variable significantly.
Percentile k1 k2 k3 k4 25 th percentile 0.3060 0.2864 0.1005 0.1111 37.5 th percentile 0.3723 0.2804 0.0831 0.0604 50 th percentile 0.4514 0.2844 0.0770 0.0608 62.5 th percentile 0.5103 0.3160 0.0721 0.0660 75 th percentile 0.5757 0.3877 0.0767 0.0744 Mean of objectives 0.5097 0.3150 0.0730 0.0653 Table III: MPSA results: ranks of all parameters for variable x at various percentile objectives. Variables k1 k2 k3 k4 x 0.0942 0.3389 0.0235 0.0870 y 0.5097 0.3150 0.0730 0.0653 z 0.6505 0.0211 0.0205 0.0958 u 0.1996 0.2837 0.1135 0.3681 Table IV: MPSA results: mean of the objectives for all parameters with variables.
From the previous studies [3,4], the melting of HbS polymers into monomers is beneficial, if not mandatory to eliminate the sickling of the red blood cell. Furthermore, carbon monoxide (CO) has been shown to improve the process of HbS polymers melting. Sangart Inc. (San Diego, CA) [19] has developed a drug called MP4CO that delivers CO at a therapeutic level to the sickle cells to facilitate the melting process. In the extended model [4], we see that many parameters play important roles in the breaking down of the de-oxy HbS polymers which leads to the formation of the CObound HbS monomers. Thus, it is crucial to determine the most important parameters that affect de-oxy HbS polymers and the CO-bound HbS monomers. Therefore, we analyze the sensitivity of all parameters to identify the most important parameters in the de-polymerization process.
Our first set of numerical experiments addresses the sensitivity of the variables with respect to the parameters using the TSA method. The representatives of sensitivity graphs are plotted as functions of time in Fig. 6 for CO monomers. The sensitivity values using the L 2 norms of the TSFs are displayed in Table I, which shows the strong effect of the CObinding rate to de-oxygenated polymers (k p (CO) on all the variables, i.e., the most sensitive parameter. The melting rate of CO-bound polymers (k CO ) shows a weak effect on the variables, thus it appears to be the least sensitive or insensitive. Therefore, the concentration of de-oxygenated HbS polymers and the concentration of CO-bound HbS monomers are affected the most by the rate k p (CO) and the least by k CO .
In the next analysis, using the previous TSFs, we study the sensitivities while all parameters are varied. As all the parameters are varied simultaneously, we notice quite a change in the sensitivity rankings. Here, we notice the importance of melting rate of de-oxygenated polymers (k d ) on the concentration of de-oxygenated HbS polymers and CO-bound HbS monomers. It should be noted that the k d is the most important parameter in the breaking down of the de-oxygenated HbS polymers and the formation of the CO-bound HbS monomers. Followed by k d , k p (CO) is the second sensitive parameter to the concentration of de-oxygenated HbS polymers and k m (CO) is the second sensitive parameter to the concentration of CO-bound HbS monomers. As for the concentration of de-oxy HbS monomers, k m (CO) followed by k d are the two most sensitive parameters. For the concentration of CO-bound HbS polymers k p (CO) is the parameter that causes the most disturbance. It should be noted that k CO has demonstrated the weakest effect in the sensitivity of all output variables.
Lastly, we perform the MPSA on the CO-mediated sickle cell de-polymerization with the same set of initial conditions and parameter ranges. The sensitivity rankings of all output values with respect to all input parameters are obtained by the MPSA analysis directly. These results are essentially identical to the results from TSA with semi-normalization of the TSFs.
In comparing the methods of sensitivity analysis employed in this study, both TSA and MPSA have pros and cons. The TSA method is simple to derive mathematically and simple to implement in the program codes to obtain the TSFs numerically. But it is limited by its large computational cost. In a case where there are many parameters (r) and output variables (n), there will be many equations (n + n × r) for the TSFs. The TSA also focuses on the local effect of the parameter on the output variables. As for the MPSA method, multiple parameters can be considered at the same time. This method studies the overall effect of the parameters on the output's variables. The disadvantage of this method is the large size of sampling for the parameters because of the Monte-Carlo simulation.
In chemical reactions, temperature plays an important role. The fluctuation of temperature influences the rate of reactions. The rate constants present in the model are influenced by the temperature variation [20]. Sensitivity analyses, such as the TSA and MPSA performed in this study, are useful tools in the iterative mathematical modeling process. Results from the sensitivity analysis can be used to inform and design future experiments, whose results can be used to further refine the mathematical model. In this study for CO-mediated sickle cell de-polymerization, the results indicate that the system shows the most sensitivity within the first 200 minutes. The system seems to approach the same equilibria regardless of the parameters chosen, which suggests that additional experimentation should focus on obtaining data during the first 200 minutes. Furthermore, this study suggests that the system is most sensitive to the melting rate of the de-oxygenated polymer (k d ),
Table V: Ranking results of TSA and MPSA.
followed by the CO binding rates for monomers and polymers (k m (CO), k p (CO)). The system shows little sensitivity to the melting rate of the CO polymer (k CO ). Future experimentation that investigates these rates would be useful to further verify these results.
For the future topics, first of all, both the TSA and MPSA analyses can be conducted for the initial conditions C d m (0), C CO m (0), C d p (0), and C CO p (0), and the solubilities C s , C CO s , for each output variable C d m (t), C CO m (t), C d p (t), and C CO p (t). Secondly, it can be seen from the analysis in this study that the sensitivity values calculated by the L 2 norm of the TSFs are affected by the scales of the functions. Instead of the TSFs, we may use the relative sensitivity functions for a more delicate study. More details of the relative sensitivity functions can be found in [5,6]. Thirdly and lastly, the MPSA can be further improved. For example, the K-S measuring is a rough approximation for the distance between two CDFs. A more delicate technique may be developed to measure the difference between the acceptance and non-acceptance distributions. Another example, if there exists appreciable correlation between parameters, the current version MPSA may not be efficient. Taking the correlation of parameters into consideration, we may project the distributions onto new axes to obtain a more accurate result.
The analysis in this study has been done in Matlab. The whole set of Matlab codes are uploaded and published on GitHub at the following link: https://github. com/lipingliuncat/ParameterSensitivity.
This study is supported by the Department of Mathematics and Statistics at North Carolina Agricultural and Technical State University. The authors greatly appreciate the detailed comments and suggestions made by the editor and reviewers, which significantly improved the work. VIII. APPENDIX The sensitivity equations for all four parameters: dx(t) dt = -k 1 x(t) C 1 -1 z(t) -k 2 x(t) d dt ∂x ∂k 1 = 1 -x(t) C 1 z(t) + k 1 1 -x(t) C 1 ∂z(t) ∂k 1 -k 1 1 C 1 ∂x(t) ∂k 1 z(t) -k 2 ∂x(t) ∂k 1 d dt ∂x ∂k 2 = k 1 1 -x(t) C 1 ∂z(t) ∂k 2 -k 1 1 C 1 ∂x(t) ∂k 2 z(t) -k 2 ∂x(t) ∂k 2 -x(t) d dt ∂x ∂k 3 = k 1 1 -x(t) C 1 ∂z(t) ∂k 3 -k 1 1 C 1 ∂x(t) ∂k 3 z(t) -k 2 ∂x(t) ∂k 3 d dt ∂x ∂k 4 = k 1 1 -x(t) C 1 ∂z(t) ∂k 4 -k 1 1 C 1 ∂x(t) ∂k 4 z(t) -k 2 ∂x(t) ∂k 4 dy(t) dt = k 3 1 -y(t) C 2 u(t) + k 2 x(t) d dt ∂y ∂k 1 = k 3 1 -y(t) C 2 ∂u(t) ∂k 1 -k 3 1 C 2 ∂y(t) ∂k 1 u(t) + k 2 ∂x(t) ∂k 1 d dt ∂y ∂k 2 = k 3 1 -y(t) C 2 ∂u(t) ∂k 2 -k 3 1 C 2 ∂y(t) ∂k 2 u(t) + k 2 ∂x(t) ∂k 2 + x(t) d dt ∂y ∂k 3 = 1 -y(t) C 2 u(t) + k 3 1 -y(t) C 2 ∂u(t) ∂k 3 -k 3 1 C 2 ∂y(t) ∂k 3 u(t) + k 2 ∂x(t) ∂k 3 d dt ∂y ∂k 4 = k 3 1 -y(t) C 2 ∂u(t) ∂k 4 -k 3 1 C 2 ∂y(t) ∂k 4 u(t) + k 2 ∂x(t) ∂k 4 dz(t) dt = k 1 x(t) C 1 -1 z(t) -k 4 z(t) d dt ∂z ∂k 1 = x(t) C 1 -1 z(t) + k 1 x(t) C 1 -1 ∂z(t) ∂k 1 + k 1 1 C 1 ∂x(t) ∂k 1 z(t) -k 4 ∂z(t) ∂k 1 d dt ∂z ∂k 2 = k 1 x(t) C 1 -1 ∂z(t) ∂k 2 + k 1 1 C 1 ∂x(t) ∂k 2 z(t) -k 4 ∂z(t) ∂k 2 d dt ∂z ∂k 3 = k 1 x(t) C 1 -1 ∂z(t) ∂k 3 + k 1 1 C 1 ∂x(t) ∂k 3 z(t) -k 4 ∂z(t) ∂k 3 d dt ∂z ∂k 4 = k 1 x(t) C 1 -1 ∂z(t) ∂k 4 + k 1 1 C 1 ∂x(t) ∂k 4 z(t) -k 4 ∂z(t) ∂k 4 -z(t) du(t) dt = -k 3 1 -y(t) C 2 u(t) + k 4 z(t) d dt ∂u ∂k 1 = -k 3 1 -y(t) C 2 ∂u(t) ∂k 1 + k 3 1 C 2 ∂y(t) ∂k 1 u(t) + k 4 ∂z(t) ∂k 1 d dt ∂u ∂k 2 = -k 3 1 -y(t) C 2 ∂u(t) ∂k 2 + k 3 1 C 2 ∂y(t) ∂k 2 u(t) + k 4 ∂z(t) ∂k 2 d dt ∂u ∂k 3 = -1 -y(t) C 2 u(t) -k 3 1 -y(t) C 2 ∂u(t) ∂k 3 + k 3 1 C 2 ∂y(t) ∂k 3 u(t) + k 4 ∂z(t) ∂k 3 d dt ∂u ∂k 4 = -k 3 1 -y(t) C 2 ∂u(t) ∂k 4 + k 3 1 C 2 ∂y(t) ∂k 4 u(t) + k 4 ∂z(t) ∂k 4 + z(t)
Biomath 13 (2024), 2312036, https://doi.org/10.55630/j.biomath.2023.12.036