Research article  Open  Published:
A qualitatively validated mathematicalcomputational model of the immune response to the yellow fever vaccine
BMC Immunologyvolume 19, Article number: 15 (2018)
Abstract
Background
Although a safe and effective yellow fever vaccine was developed more than 80 years ago, several issues regarding its use remain unclear. For example, what is the minimum dose that can provide immunity against the disease? A useful tool that can help researchers answer this and other related questions is a computational simulator that implements a mathematical model describing the human immune response to vaccination against yellow fever.
Methods
This work uses a system of ten ordinary differential equations to represent a few important populations in the response process generated by the body after vaccination. The main populations include viruses, APCs, CD8+ T cells, shortlived and longlived plasma cells, B cells and antibodies.
Results
In order to qualitatively validate our model, four experiments were carried out, and their computational results were compared to experimental data obtained from the literature. The four experiments were: a) simulation of a scenario in which an individual was vaccinated against yellow fever for the first time; b) simulation of a booster dose ten years after the first dose; c) simulation of the immune response to the yellow fever vaccine in individuals with different levels of naïve CD8+ T cells; and d) simulation of the immune response to distinct doses of the yellow fever vaccine.
Conclusions
This work shows that the simulator was able to qualitatively reproduce some of the experimental results reported in the literature, such as the amount of antibodies and viremia throughout time, as well as to reproduce other behaviors of the immune response reported in the literature, such as those that occur after a booster dose of the vaccine.
Background
Mathematical and computational modeling is constantly evolving tool, which can be applied to many distinct research areas, such as Biology, Physics, Chemistry, Engineering, Biomechanics, Climate Modeling, tsunami and earthquake prediction, among others [1–23]. With this type of tool, the phenomenon under study is represented by mathematical equations which can be solved using computational simulators. The use of such mathematicalcomputational models can help reduce costs, time, risks and volunteers involved in the research. However, to achieve these objectives, the models must be very reliable. Since models are always an abstraction of reality, using simplifications to deal with complexities, many factors that can contribute to the real phenomenon may be ignored.
Mathematical models have been used for many years to represent various aspects of the immune system and related pathologies, but their application to describe the effects of vaccines has been rather limited [24]. The term computational vaccinology has been used to refer to computeraided vaccine design [25–28], and its objective is to use different modeling techniques to aid the development and improvement of vaccines at different stages of their design processes.
In 1796 Edward Jenner introduced vaccination against smallpox, which was a major health problem at the time. Jenner observed that milkmaids were protected from smallpox after having suffered from cowpox, and concluded that cowpox could be used as a deliberate mechanism of protection against smallpox [29]. Jenner inoculated an 8yearold boy, James Phipps, with cowpox. Subsequently, Jenner inoculated the boy again, this time with smallpox, and he did not contract the disease [29]. Jenner concluded that protection was complete. This is the key of vaccination: expose the body to antigens from pathogens, in order to stimulate the production of antibodies and defense cells against a specific disease.
Much of the work on computational vaccinology is related to the process of creating a new vaccine, such as in the selection of the best strains for use. In a previous paper [30] we proposed a new use for computational vaccinology, i.e. in the clinical development stage. With the use of mathematical and computational models, it is possible to experiment, in silico, different scenarios related to vaccination, to address important questions that remain unanswered.
The Yellow Fever (YF) vaccine, available since 1937 [31], is made from live attenuated virus still capable of triggering an immune response and inducing the production of antibodies and memory cells. Livevirus vaccines induce an immune response similar to that obtained with exposure to wild virus, but the risk of presenting characteristic symptoms of the disease and its complications, or death, due to vaccination, is extremely small. The YF vaccine is considered an effective and safe vaccine, with high documented seroconversion rates and low rates of adverse events. It has been effectively used to control a noneradicable disease and is one of the vaccines that can benefit from the new use of computational immunology. The reason for this is manifold.
Although recognized as an effective and safe vaccine, some questions remain unanswered or poorly understood, and could be reassessed using new technologies and tools. As the vaccine was developed decades ago, some steps of its developmental processes were established empirically. A good example is the optimal dose required for immunization. What is the vaccine dose with the best immunogenicity/reatogenicity ratio? There are clinical studies designed to evaluate this [32], but these studies require time and resources, and there are methodological restrains to test several different doses. With the use of mathematical and computational modeling techniques, it is possible to evaluate a larger spectrum of doses in a much shorter time, using far less resources.
Another controversial issue is the need for a booster dose. Using mathematical and computational modeling, it is possible to simulate, for an individual, what his/her antigen levels will be years in the future, in a few minutes, to assess the duration of immunity and the need for booster dose administration, taking into account differences among individuals and doses,to help in the design of prospective cohort studies.
Despite being considered a safe vaccine, there are rare serious adverse events that need to be reassessed, such as viscerotropic and neurotropic events. There are also questions regarding the safety for vaccinating specific populations such as the elderly, people living with Human Immunodeficiency Virus (HIV)/AIDS and other immunocompromised populations. Because the YF vaccine is a liveattenuated virus vaccine, there is a small but not insignificant risk of occasional higher viral replication related either to vaccine virus attenuation aspects or an inability of the immune system to control the vaccine virus replication.
Recently, YF outbreaks were recorded in Angola and the Democratic Republic of Congo (DRC), with the latest outbreak still underway in Brazil, starting in December 2016. From December 2016 to February 22, 2017, 1,345 suspected cases were recorded, of which 295 have been confirmed, and 215 deaths reported to the Brazilian Ministry of Health [33].
YF is not an eradicable disease because of its sylvatic cycle. Reported cases in the Brazilian outbreaks were classified as sylvatic YF, but the risk of urban YF reintroduction is imminent due to high levels of Aedes aegypti infestation in Brazilian cities where vaccination coverage is not routinely recommended. World stocks and YF vaccine production capacity are a logistic concerns which could impact the control of disease transmission, particularly in large outbreaks. In Kinshasa, capital of the DRC, fractional doses of the YF vaccine were administered for outbreak control [34].
Concerned about the risk of a global epidemic, the WHO launched in April 2017 a strategy called Eliminate Yellow fever Epidemics (EYE), which aims to eliminate YF epidemics in the world by 2026 [34]. Through early detection and rapid and appropriate response, it is possible to minimize suffering, damage and propagation [34]. This strategy has three goals: protect populations at risk, prevent the international spread of YF and contain outbreaks quickly. To achieve these goals, the strategy suggests actions on different fronts, including research and development of better tools and practices. Assessing data about optimal vaccine dose and duration of immunity could help the design of new vaccination strategies for disease control.
This work, therefore, presents a first step towards an ideal scenario to simulate distinct situations related to the use of the YF vaccine: a qualitatively validated mathematicalcomputational model of the immune response to the YF vaccine. The model considers the major populations of Human Immune System (HIS) cells and molecules important in the process of immunity acquisition, such as Antigen Presenting Cells (APCs), B and T lymphocytes, and antibodies, which are considered the main marker of immunity. The model was then evaluated using distinct scenarios, and was successful in qualitatively reproducing experimental results reported in the literature.
This work is organized as follows. First, Section Related works presents related works done in this fields. Section Methods presents the mathematical and computational models used to reproduce the immune system response to the YF vaccine. The results are then presented in Section Results and discussed in Section Discussion, and finally Section Conclusion presents our conclusions.
Related works
The use of mathematical and computational models to help vaccine development is not new. In fact, several works use computational tools to aid vaccine design. For example, epitopemapping algorithms have been used for vaccine design since the 1980s [35]. Since then, new computational tools have been used for selection of vaccine targets [36–44]. Most of the works focuses on using mathematical and computational tools to predict epitopes [45] or to develop virtual screening approaches (i.e, the identification of relevant antigens) [46–49]. This traditional use of computational vaccinology is related to preclinical development. This work focuses on the development of mathematical and computational models that can be used in the clinical development stage, i.e., when the vaccine is first tested in humans. We argue that it is possible to carry out some experiments in silico, reducing the search space for experiments in vivo or in vitro, and it is possible to eliminate, reinforce or weaken hypotheses and to propose new studies, thus saving time and resources.
Several computational modeling techniques applied to vaccination are analyzed and discussed by Pappalardo et al. [24]. The authors describe what mathematical and computational modeling are and how they can aid research in vaccination. Modeling is defined as human activity involving the representation, manipulation, and communication of everyday real world objects. In their review, two main types of modeling are considered: AgentBased Models (ABM) and mathematical models. Mathematical models are mainly based on differential equations, whether ordinary or partial, with delayed and/or stochastic equations. In this work, we try to qualitatively validate a simplified mathematicalcomputational model of the immune response to the YF vaccine presented in a previous work [50], which is based on a live, attenuated viral strain. The model uses Ordinary Differential Equations (ODEs) to model the main cells and molecules related to adaptive immune response.
Another work uses an ODEbased approach to model the human immune response to vaccination against both YF and smallpox [51] using distinct data and equations sets, one for each disease. The aim of the authors was to primarily evaluate the dynamics of CD8+ T cells, while this work will evaluate the immune response as a whole. The model proposed here differs from that presented by Le et al. [51], since it considers important populations at each stage of the immune response to YF vaccination, from virus inoculation to APC antigen presentation and consequent activation of lymphocytes, generation of antibodies and memory cells.
Methods
Mathematical model
In this section, we present the model we proposed in a previous work [50], which will be qualitatively validated in this paper. The model consists of a system of 10 ODEs representing important populations in the response process generated by the body after vaccination. The main populations are viruses, APCs, CD8+ T cells, shortlived and longlived plasma cells, B cells and antibodies. Only populations related to the YF vaccine are modeled. For example, only B and T cells whose receptors can recognize the YF virus are considered in the model.
Equation 1 represents the vaccine virus (V).
The virus can not proliferate by itself, it needs to infect a cell and use it as a factory for new viruses. This is implicitly considered in the term π_{ v }V, which represents the multiplication of the virus in the body, with a production rate of π_{ v }. The term $\frac {c_{v1} V}{c_{v2} + V} $ denotes a nonspecific viral clearance by the innate immune system. This function is similar to the Hill family of equations [52].
The above equation is a generalization of the hyperbolic saturation function. The parameter k_{1} scales the maximum value to which the function is asymptotic, k_{2} is a shape parameter and k_{3} is analogous to the halfsaturation constant. If k_{2}=1, the MichaelisMenten function is produced [53].
The term k_{v1}VA denotes specific viral clearance due to antibody signaling, where k_{v1} is the clearance rate. The term k_{v2}VT_{ E } denotes specific viral clearance due to the induction of apoptosis of cells infected by the YF virus, where k_{v2} is the clearance rate.
APCs are all cells that display antigens complexes on their surfaces, such as dendritic cells and macrophages. Two stages of APCs were considered: immature and mature. The first stage, immature APCs (A_{ P }), is described by Eq. 2.
The term $\alpha _{A_{P}} \left (A_{P0}  A_{P}\right)$ denotes the homeostasis of APCs, where $\alpha _{A_{P}}$ is the homeostasis rate. The term $\beta _{A_{P}} A_{P} \left (k_{A_{P1}} + tanh\left (V  k_{A_{P2}}\right)\right) $ denotes the conversion of immature APCs into mature ones. Therefore, the same term appears in Eq. 3 with positive sign. The constant $\beta _{A_{P}}$ represents the conversion rate and $\left (k_{A_{P1}} + tanh\left (V  k_{A_{P2}}\right)\right)$ is a sigmoidal saturation function in the form of a hyperbolic tangent.
Equation 3 represents the mature APCs (A_{ PM }).
The first term, as explained, denotes the dynamics of APCs maturation. The second term, $\delta _{A_{PM}}A_{PM}$, denotes the natural decay of the mature APCs, where $\delta _{A_{PM}}$ is the decay rate.
Equation 4 represents the population of naïve CD8+ T cells (T_{ N }).
The term $\alpha _{T_{N}} (T_{N0}  T_{N})$ represents the homeostasis of CD8+ T cells, where $\alpha _{T_{N}}$ is the homeostasis rate. The term π_{ T }A_{ PM }T_{ N } denotes the activation of naïve the CD8+ T cells, where π_{ T } is the activation rate. Therefore, the same term appears in Eq. 5 with positive sign.
Equation 5 represents the effector CD8+ T cell population (T_{ E }).
The term $k_{T_{E}} A_{PM} T_{E}$ represents the proliferation of effector CD8+ T cells. The term $\delta _{T_{E}} T_{E} $ represents the natural death of these cells, with $\delta _{T_{E}}$ representing its decay rate.
Equation 6 represents B cells (B), both naïve and effector ones. These populations were not considered separately in order to simplify the model.
The term α_{ B }(B_{0}−B) represents the B cells homeostasis, where α_{ B } is the homeostasis rate. The term π_{ B }A_{ PM } represents the proliferation of the active B cells. The terms β_{ S }A_{ PM }B, β_{ L }A_{ PM }B and $\beta _{B_{M}} A_{PM} B$ denote the portions of active B cells that differentiate into shortlived plasma cells, longlived plasma cells and memory B cells, respectively. These terms will appear with positive sign in Eqs. (7), (8) and (9). The activation rates are respectively given by β_{ S }, β_{ L } and $\beta _{B_{M}}$.
Equation 7 represents the shortlived plasma cells (P_{ S }).
The term δ_{ S }P_{ S } denotes the natural decay of shortlived plasma cells, where δ_{ S } is the decay rate.
Equation 8 represents the longlived plasma cells (P_{ L }).
The term δ_{ L }P_{ L } denotes the natural decay of longlived plasma cells, with δ_{ L } representing the decay rate. The term γ_{ M }B_{ M } represents the production of these cells by memory B cells, where γ_{ M } is the production rate.
Equation 9 corresponds to memory B cells (B_{ M }).
The term $k_{B_{M1}} B_{M}\left (1  \frac {B_{M}}{k_{B_{M2}}}\right)$ represents the logistic growth of memory B cells, i.e., there is a limit to this growth. The constants $k_{B_{M1}}$ and $k_{B_{M2}}$ represent the growth rate and limits, respectively.
Equation 10 represents the antibodies. The terms π_{ AS }P_{ S } and π_{ AL }P_{ L } are the production of the antibodies by shortlived and longlived plasma cells, respectively. The production rates are given by π_{ AS } and π_{ AL }, respectively. The term δ_{ A }A denotes the natural decay of these cells, where δ_{ A } is the decay rate.
The model presented in this paper was based on an earlier study [54], which described a mathematical model to represent the human immune response to an infection by YF virus. Therefore, the first difference is that this paper focus on modelling the effects of the YF vaccine administered subcutaneously.
The previous work [54] modeled the immune response to the YF virus from infection of epithelial cells to secretion of antibodies, considering various populations of cells and molecules, in different stages and compartments. There were 19 ODEs divided into two compartments: one representing the tissue where the virus proliferates and the other the lymph nodes. In order to consider all the cells and molecules, the model became complex.
Another issue is related to its adjustment to reproduce some behaviors described in the literature: as the number of equations and parameters increases, so does the amount of data and information needed to adjust the model. The second difference between the two models is that the model reproduced here [50] reduces the number of equations from 19 to 10. The reduced model reproduced in this work considers only the main populations of cells and molecules involved in the response to the vaccine, and abstracts some details that are not crucial to represent the behavior of the immune response, such as the representation of distinct compartments. In addition, some populations were not considered because no experimental data are available to validate the simulations, such as CD4+ T cells. In the near future, more cells or molecules can be reintroduced in the model if their roles are important to explain or represent behaviors that the reduced model [50] could not represent. Table 1 summarizes the main differences between the model presented in previous work [54] and the one evaluated in this work.
It is important to remember that a mathematical model is an abstraction of reality and therefore simplifications are always necessary. This is accentuated when the target of the model is the HIS response, a complex network that involves many tissues, organs and cells and that performs several processes. The level of abstraction depends on the purpose of the model. HIS can be seen at various levels, from the level of substances produced by cells, such as cytokines, to the level of cells and molecules, as in the case of the simplified model [50]. It also can reach the level of an entire population, as in the case of the epidemiological models. The use of a simplified model does not imply that it can reproduce only a limited number of scenarios. The point is that some of the aspects not directly included in the model may be indirectly present, as constants. As such, the choice of distinct values for some constants may represent distinct behaviors in the system.
Computational model
For the resolution of the ODEs system, a code was implemented using Python programming language, which includes libraries for easily solving complex mathematical problems. The library chosen was SciPy [55]. This library has a package called “integrate”. One of the functions available in this package is called “odeint”, and it is used to numerically solve a system of ODEs. The choice of the numerical method to be used is made automatically by the function based on the characteristics of the equations. The function uses an adaptive scheme for both the integration step and the convergence order. The function can solve the ODEs system using either the Backward Differentiation Formula (BDF) or the Adams method [56]. BDF is used for stiff equations and the implicit Adams method is used otherwise.
The experiments were performed using Python version 2.7.10 using the Spyder Integrated Development Environment (IDE). The execution environment was composed by an Intel Core i5 1.6 GHz processor, with 8 GB of RAM. The system runs macOS Sierra version 10.12.5.
Results
In order to qualitatively validate our model, four experiments were carried out. The first one simulates a scenario where an individual was vaccinated against YF for the first time. The standard dose of the vaccine was used in this scenario. The results of the simulation were then compared to experimental data obtained from the literature.
The second scenario assesses the immune response following the administration of a booster dose ten years after the first dose. The standard dose of the vaccine was also used in this scenario. Although there are no experimental data from the literature that could be used for comparison purposes, there is research reporting that the expected behavior is an individual to present a lower viremia, and to raise antibodies levels to levels higher than those obtained after the administration of the first dose [57].
The third scenario simulates the immune response to the YF vaccine in individuals with different levels of naïve CD8+ T cells prior to vaccination. This simulation aims to evaluate the importance of this population of cells in the control of viremia and in the production of antibodies.
The fourth scenario is based on an experimental study [32] in which distinct doses of the YF vaccine were tested. Compared to the standard dose, the experimental study reported that, to some extent, the reduction did not significantly affect the percentage of seroconversion. In this scenario, computational experiments are executed several times, using distinct values for the vaccine doses. For comparison purposes, the computational experiments were carried out using the same values of the experimental study [32].
In addition to evaluating the response of the model when different doses are administered, we performed a sensitivity analysis of the parameters related to virus dynamics. The main results of this analysis are shown in subsection Sensitivity analysis.
All the initial values used for the variables as well as the model parameters are presented in Tables 2 and 3. The parameters were adjusted, except for δ_{ A }, whose value was extracted from the literature.
In general, the literature reports two distinct sets of experimental data. The first one is viremia along time, i.e., the amount of virus present in the bloodstream. The second dataset reported in the literature is the antibody levels along time. Therefore, in order to validate the model, the values obtained by Eqs. 1 and (10) will be compared to experimental values found in the literature.
First vaccination
This section presents the computational results of a simulation in which an individual was vaccinated against YF for the first time.
In this computational experiment, a value equal to 27,476 International Units (IU) was used as the standard amount of virus particles present in the vaccine. This value is set as the initial condition of the virus population represented in Eq. 1 by V (all initial conditions are presented in Table 2). In fact, in the case of the 17DDYFV the amount of virus particles varies depending on the vaccine lot number, ranging from 2.3 to 12 times the minimum value required by the WHO [32]. The 17DDYFV is the YF vaccine developed by BioManguinhos/Fiocruz, one of the three producers prequalified by the WHO to supply vaccines to international agencies.
Figures 1 and 2 show the comparison of the antibody curve generated as a result at 100 and 4,000 simulation days, respectively, with the experimental results from the literature [58]. The result of the simulation is presented in separate figures in order to better observe the increase of the antibody level in the first days after vaccination.
The levels of antibodies obtained from the literature [58] are in Geometric Mean Titers (GMT) and refer to time intervals after vaccination. The time values used in the graph were obtained by averaging the times of each interval. For example, the first point was the 3045 days postvaccination interval, the value used was 37 days, the corresponding antibody level was 8,762.8 IU/mL.
Figure 3 shows the viremia curve obtained by the simulation of the model in comparison with the experimental data obtained from the literature [32].
Booster dose
The administration of a booster dose was simulated 10 years after the administration of the first dose. The simulation is quite simple. As in the previous scenario, the initial value of V was set to 27,476 to simulate the administration of the first dose. Then, the simulation is executed until day number 3,650, when the value of variable V is set again to 27,476. The difference from the beginning of the simulation is that this time antibodies and memory cells that were produced after the first dose are present. Figures 4 and 5 present the antibody curves. The results of the booster dose simulation were shifted to facilitate its comparison to the results of the first dose. The solid line curve represents the response to the first dose, while the dashed curve represents the response to the booster dose.
Figure 6 shows the viremia curves 15 days after the administration of the vaccine.
Naïve T CD8+ levels
The clearing of the intracellular pathogen via CD8+ cytotoxic T lymphocytes appears to be important for recovering from primary viral infection. Based on this observation found in the literature [57], the authors decided to compare the immunological response given by the simulation of the model with different levels of CD8+ T cells in order to evaluate the impact of this population of lymphocytes on viral clearance.
Figure 7 shows the viremia curves for different initial CD8+ T cell values, and Figs. 8 and 9 show antibody levels.
Doseresponse
An experimental work [32] reported that “doses from 27,476 IU to 587 IU induced similar seroconversion rates and neutralizing antibodies geometric mean titers (GMTs)”. Based on this study, a second scenario analyzes the results of our model when different dose values are administered. The values used in the simulation are the same to those used by the experimental work [32]: 31 IU, 158 IU, 587 IU, 3,013 IU, 20,447 IU and 27,476 IU.
Figures 10 and 11 show the viremia curves obtained by the model for distinct vaccine doses. Figure 11 uses a smaller scale to allow the visualization of the simulated viremia curve obtained after administration of the dose with 587 IU (represented by diamonds).
Figure 12 presents the antibody curves generated by the computational model during a 50day period for different doses of the vaccine, while Fig. 13 presents the antibody curves obtained by simulating 4,000 days after vaccination.
Doses using 31 IU and 158 IU did not produce viremia nor antibody titers, so the curves are superimposed on the xaxis.
Sensitivity analysis
The sensitivity analysis identifies the impact caused by the variation of parameters and initial conditions of the mathematical model in the dependent variables [59]. If a small change in a parameter is responsible for a drastic change in the result of the problem, it means that the problem is sensitive to that particular parameter. Otherwise, this parameter has a low impact on the model. This analysis is used to help the understanding of the mathematical model, since it allows the identification of the most relevant parameters, that is, the values of these parameters must be carefully defined.
A bruteforce approach was used to examine the influence of all parameters of the model. The parameter values were varied from 10% to + 10% (in 5% intervals) from their original values. The original values are presented in Table 3 and were obtained after adjustment using experimental data from the literature [32, 58, 60, 61]. For each parameter, the curves that simulate the level of antibodies and viremia were evaluated, since they are the main populations of interest and on which there is experimental data. Only the parameters to which the model was most sensitive will be presented in this section.
As expected, the model was sensitive to most of the parameters of the equation describing the virus dynamics (Eq. 1). The antibody curves were not significantly affected by them, therefore only the viremia curves will be presented. Figure 14 shows the distinct viremia curves obtained for different values of π_{ v }. This parameter represents the viral replication rate and, as one could expect, the model was very sensitive to it. The more the virus can multiply, the more difficult it is for the HIS to contain it and the higher the viremia level is.
Figures 15 and 16 present the viremia curves obtained by simulation of the model for different values of k_{v1} and k_{v2}, respectively. These parameters represent the neutralization rates of the YF virus per unit of neutralizing antibodies (k_{v1}) and CD8+ T cells (k_{v2}).
If we consider that the parameter k_{v1} represents the ability of antibodies to neutralize the YF virus, its value can be understood as its affinity/specificity to the YF virus and, if it is more specific and can neutralize the virus better, viral replication will be better controlled and viremia will be lower.
It is easy to understand why the model, especially the viremia curve, is so sensitive to parameter k_{v2}. It represents the ability of CD8+ T cells to induce apoptosis of an infected cell. Thus, the higher this ability, the fewer the number of infected cells. Since YF viruses use infected cells to reproduce themselves, the viremia level is reduced.
Discussion
As can be observed from Figs. 1 and 2, from a qualitative point of view, the values obtained from the computational experiments are very close to the experimental results. Also, the literature reports that the antibody concentration in the bloodstream peaks at about two weeks after vaccination [62], a value close to the one obtained in the computational experiments.
Figure 3 shows that, in the simulation, the peak viremia value occurs on the fifth day, consistent with the literature, which reports that it occurs between four and six days after vaccination [63], as well as with experimental results [32]. In addition, the literature reports that ten days after vaccination, viremia is undetectable [63], which is consistent with the computational results. For some patients, however, viremia can be detectable, as one experimental result has shown [32].
Figures 4 and 5 show that the behavior described in the literature [57] resembles that obtained by the simulation of the model: the neutralizing antibodies levels are slightly increased after the booster dose.
The viremia curves shown in Fig. 6 demonstrate that viremia reaches much lower levels after the administration of the booster dose than those seen after the administration of the first dose. Although the viremia level is lower after the administration of the booster dose than the first dose, it is not possible to say that this level is below the threshold of detectable viremia as described in the literature: “viremia has not been documented in persons receiving a booster dose of YF vaccine” [64]. This occurs because of the use of distinct units to measure viremia, and the fact that it is not trivial to convert one unit to another. For this reason, this work considers only qualitative results, and not quantitative ones. This result still needs to be quantitatively compared to experimental data in order to better validate our model, but the qualitative behavior presented is satisfactory since the level of antibodies and/or memory cells was able to contain viral replication more efficiently than it was observed for a naïve individual. Furthermore, the viremia level was almost 4 times lower for the booster dose than for the first dose.
As show in Fig. 7, as the number of CD8+ T cells is reduced, the viremia increases and lasts longer, reinforcing the importance of CD8+ T cells in the control of viral replication. Figures 8 and 9 show that this variation in initial CD8+ T cell values did not significantly affect antibody production nor duration of immunity.
As observed in the Figs. 10 and 11, all doses greater than 3,013 IU produce high levels of viremia. Although the viremia increases with the use of doses with higher concentrations, the antibody response presents a very subtle difference, as can be observed in Fig. 12. The 587 IU dose, which presented a much smaller, unremarkable viremia (Fig. 10) compared to the doses with higher concentrations, was also able to induce an antibody response similar to that induced by formulations with higher concentrations.
Figure 13 suggests that the duration of immunity does not appear to be affected by vaccine formulations with distinct concentrations: all doses above 587 IU present similar results. For now, it seems that yellow fever vaccine can be used in much lower doses than usual: the computational experiments indicate that vaccine formulations with 587 IU can produce the same seroconversion rates than the 27,476 IU formulation, in accordance to the experimental data [32]. Although the reference paper investigated the duration of immunity for a smaller period of time, approximately 10 months after vaccination [32], its conclusions were similar to those obtained by the computational experiments: “GMTs of each group were not statistically different from the reference vaccine". The computational results are also in agreement with other studies. One paper [57] concluded: “there was no correlation between the level and duration of detectable 17D viremia and the postvaccination nAb level". Another paper presents a similar conclusion [65]: “the serological response was not related to virus dose as the titres obtained with high or low doses of virus was at the same level".
The results of the simulations for these four scenarios have shown that the model was able to reproduce, from a qualitative perspective, clinical results reported in the literature, despite all simplifications [50]. Some aspects not directly included in the model may be indirectly present, as constants. Therefore, the choice of distinct values for a constant may represent distinct behaviors in the system. For example, one paper [58] points out that “The decreasing trend in antibody titres with the time since vaccination appeared strongly modified by age". In our model, the effects of age in the production of antibodies could be reproduced increasing or decreasing the values used for the antibody secretion rate (Eq. 10, π_{ AS } and π_{ AL }). For this reason, our model does not need to include age as one of its parts. The same applies for other aspects of the immune system that are not directly included in the set of equations.
In this work we consider that the vaccine does not cause adverse events, such as Yellow fever vaccineassociated viscerotropic disease (YELAVD) and Yellow fever vaccineassociated neurotropic disease (YELAND), due to their rarity.
This model was developed and adjusted based on the immune system response to the YF vaccine, but it should be noted that the concept presented in the mathematical model is generic enough to represent the action of other diseases or vaccines in the HIS. For this, changes in both initial conditions and parameters values are probably needed.
Obtaining experimental data to adjust and validate the model is not a trivial task. Studies on the duration of immunity are difficult to interpret because different groups use distinct methods to evaluate seroprotection. There is no wellestablished serological value of protection in humans and cellular immunity data are very scarce. Also, as stated above, the values reported for viremia use distinct units, which cannot be converted into other units due to the different methods used to obtain such data. These factors made it difficult to obtain experimental data compatible with the standards and units used in the model presented in this work, and consequently to use more studies available in the literature to adjust and validate it.
Although the results found are qualitatively in agreement with the few experimental data found in the literature, more tests and refinement of the model may be needed to adjust it. To do so, experimental data to better validate the simulated scenarios needs to be obtained, in particular for booster dose and CD8+ T cells. With this, it would be possible to also evaluate the model from a quantitative perspective and, if necessary, to better adjust it. With more data available, the model may be improved, making it more reliable and sufficiently accurate to be used to help answer open questions about YF vaccine.
One of the next steps in this work is to reintroduce CD4+ T lymphocytes in the model. This could be important to simulate the effects of the YF vaccine in immunosuppressed individuals, such as people living with HIV. Since many of the YF endemic countries are in Africa, where the HIV infection rates are also high, the investigation of the best YF vaccination scheme for these individuals is relevant, since they have, in general, fewer CD4+ T lymphocytes, which are important to the activation of other lymphocytes and consequently to the production of antibodies. This population deserves special attention because the YF vaccine is made with live virus, therefore the risk of systemic lethal infection exists.
Conclusion
This work presented the qualitative validation of a reduced mathematicalcomputational model to represent the immune response to the YF vaccine using four distinct scenarios. The first one simulates the immune response to the administration of the standard dose of the 17DDYFV. The second one simulates the immune response to distinct doses of vaccine. The third scenario simulates the administration of a booster dose ten years after the first dose. Finally, we evaluated the impact of changing the CD8+ T cells values. The results of a sensitivity analysis of the model was also shown. Two populations, virus and antibodies, were the main focus of the simulations because more experimental data are available and qualitative behaviors are described in the literature for these populations. The results of the simulations were collected and compared to the values reported in the literature. From a qualitative point of view, the results obtained by the computational model satisfactorily reproduced the clinical results.
Abbreviations
 ABM:

Agentbased models
 APCs:

Antigen presenting cells
 BDF:

Backward differentiation formula
 DRC:

Democratic Republic of Congo
 EYE:

Eliminate yellow fever epidemics
 GMT:

Geometric mean titers
 HIS:

Human immune system
 HIV:

Human immunodeficiency virus
 IDE:

Integrated development environment
 IU:

International units
 ODEs:

Ordinary differential equations
 YF:

Yellow fever
 YELAND:

Yellow fever vaccineassociated neurotropic disease
 YELAVD:

Yellow fever vaccineassociated viscerotropic disease
References
 1
Schoeberl B, EichlerJonsson C, Gilles ED, Müller G. Computational modeling of the dynamics of the map kinase cascade activated by surface and internalized egf receptors. Nat Biotechnol. 2002; 20(4):370–5.
 2
Wiley HS, Shvartsman SY, Lauffenburger DA. Computational modeling of the egfreceptor system: a paradigm for systems biology. Trends Cell Biol. 2003; 13(1):43–50.
 3
Doddi SK, Bagchi P. Threedimensional computational modeling of multiple deformable cells flowing in microvessels. Phys Rev E. 2009; 79(4):046318.
 4
Beard DA, Schlick T. Computational modeling predicts the structure and dynamics of chromatin fiber. Structure. 2001; 9(2):105–14.
 5
Querec TD, Akondy RS, Lee EK, Cao W, Nakaya HI, Teuwen D, Pirani A, Gernert K, Deng J, Marzolf B, et al. Systems biology approach predicts immunogenicity of the yellow fever vaccine in humans. Nat Immunol. 2009; 10(1):116–25.
 6
Karr JR, Sanghvi JC, Macklin DN, Gutschow MV, Jacobs JM, Bolival B, AssadGarcia N, Glass JI, Covert MW. A wholecell computational model predicts phenotype from genotype. Cell. 2012; 150(2):389–401.
 7
Clarke S, Vvedensky DD. Origin of reflection highenergy electrondiffraction intensity oscillations during molecularbeam epitaxy: A computational modeling approach. Phys Rev Lett. 1987; 58(21):2235.
 8
Sakurai T. Computational modeling of magnetic fields in solar active regions. Space Sci Rev. 1989; 51(12):11–48.
 9
Cuitino AM, Ortiz M. Computational modelling of single crystals. Model Simul Mater Sci Eng. 1993; 1(3):225.
 10
Yanez J, Kuznetsov M. An analysis of flame instabilities for hydrogen–air mixtures based on sivashinsky equation. Phys Lett A. 2016; 380(33):2549–2560.
 11
Feldgus S, Landis CR. Largescale computational modeling of [rh (duphos)]+catalyzed hydrogenation of prochiral enamides: reaction pathways and the origin of enantioselection. J Am Chem Soc. 2000; 122(51):12714–27.
 12
Bicerano J. Computational Modeling of Polymers. New York: CRC press; 1992.
 13
Rots JG. Computational modeling of concrete fracture. PhD thesis, Technische Hogeschool Delft. 1988.
 14
Schafer B, Peköz T. Computational modeling of coldformed steel: characterizing geometric imperfections and residual stresses. J Constr Steel Res. 1998; 47(3):193–210.
 15
Roussel N, Geiker MR, Dufour F, Thrane LN, Szabo P. Computational modeling of concrete flow: general overview. Cem Concr Res. 2007; 37(9):1298–307.
 16
McHugh P, Asaro R, Shih C. Computational modeling of metal matrix composite materials. i. isothermal deformation patterns in ideal microstructures. Acta Metallurgica et Materialia. 1993; 41(5):1461–76.
 17
Porter B, Zauel R, Stockman H, Guldberg R, Fyhrie D. 3d computational modeling of media flow through scaffolds in a perfusion bioreactor. J Biomech. 2005; 38(3):543–9.
 18
Kuhl E, Maas R, Himpel G, Menzel A. Computational modeling of arterial wall growth. Biomech Model Mechanobiol. 2007; 6(5):321–31.
 19
Randall DA, Ringler TD, Heikes RP, Jones P, Baumgardner J, et al. Climate modeling with spherical geodesic grids. Comput Sci Eng. 2002; 4(5):32–41.
 20
Nefedova V, Jacob R, Foster I, Liu Z, Liu Y, Deelman E, Mehta G, Su MH, Vahi K. Automating climate science: Large ensemble simulations on the teragrid with the griphyn virtual data system. In: 2006 Second IEEE International Conference on eScience and Grid Computing (eScience’06). Washington, DC: IEEE Computer Society: 2006. p. 32–32.
 21
Bernholdt D, Bharathi S, Brown D, Chanchio K, Chen M, Chervenak A, Cinquini L, Drach B, Foster I, Fox P, et al. The earth system grid: Supporting the next generation of climate modeling research. Proc IEEE. 2005; 93(3):485–95.
 22
Das S, Aki K. Fault plane with barriers: a versatile earthquake model. J Geophys Res. 1977; 82(36):5658–70.
 23
Loomis HG. Tsunami prediction using the reciprocal property of green’s functions. Mar Geodesy. 1979; 2(1):27–39.
 24
Pappalardo F, Flower D, Russo G, Pennisi M, Motta S. Computational modelling approaches to vaccinology. Pharmacol Res. 2015; 92:40–5.
 25
Doytchinova IA, Flower DR. Quantitative approaches to computational vaccinology. Immunol Cell Biol. 2002; 80(3):270.
 26
Brusic V, Petrovsky N. Bioinformatics for characterisation of allergens, allergenicity and allergic crossreactivity. Trends Immunol. 2003; 24(5):225–8.
 27
Taylor PD, Flower DR. In: Flower D, Timmis J, (eds).Immunoinformatics and Computational Vaccinology: A Brief Introduction. Boston: Springer; 2007, pp. 23–46.
 28
Flower DR. Bioinformatics for Vaccinology. United Kingdom: John Wiley & Sons; 2008.
 29
Paul WE. Fundamental Immunology, 5th edn. Philadelphia: Wolters Kluwer/Lippincott Williams & Wilkins; 2008.
 30
Bonin CRB, Fernandes GC, dos Santos RW, Lobosco M. Mathematical modeling based on ordinary differential equations: A promising approach to vaccinology. Hum Vaccines Immunotherapeutics. 2017; 13(2):484–9.
 31
Theiler M, Smith HH. The use of yellow fever virus modified by in vitro cultivation for human immunization. J Exp Med. 1937; 65(6):787–800.
 32
Martins RM, Maia MdLS, Farias RHG, Camacho LAB, Freire MS, Galler R, Yamamura AMY, Almeida LFC, Lima SMB, Nogueira RMR, et al. 17dd yellow fever vaccine: a double blind, randomized clinical trial of immunogenicity and safety on a doseresponse study. Hum Vaccines Immunotherapeutics. 2013; 9(4):879–88.
 33
Goldani LZ. Yellow fever outbreak in brazil, 2017. Braz J Infect Dis. 2017; 21(2):123–4.
 34
WHO. Weekly epidemiological record. 2017. http://apps.who.int/iris/bitstream/10665/255040/1/WER9216.pdf?ua=1. Accessed 13 May 2018.
 35
DeLisi C, Berzofsky JA. Tcell antigenic sites tend to be amphipathic structures. Proc Natl Acad Sci. 1985; 82(20):7048–52.
 36
Kumar N, Hendriks BS, Janes KA, de Graaf D, Lauffenburger DA. Applying computational modeling to drug discovery and development. Drug Discov Today. 2006; 11(1718):806–11.
 37
De Groot AS, Moise L, McMurry JA, Martin W. In: Falus A, (ed).EpitopeBased ImmunomeDerived Vaccines: A Strategy for Improved Design and Safety. New York: Springer; 2009, pp. 39–69.
 38
Oliveira FM, Coelho IE, Lopes MD, Taranto AG, Junior MC, Santos LL, Villar JA, Fonseca CT, Lopes DD. The use of reverse vaccinology and molecular modeling associated with cell proliferation stimulation approach to select promiscuous epitopes from schistosoma mansoni. Appl Biochem Biotechnol. 2016; 179(6):1023–40.
 39
Rappuoli R, Bottomley MJ, D’Oro U, Finco O, De Gregorio E. Reverse vaccinology 2.0: Human immunology instructs vaccine antigen design. J Exp Med. 2016; 213(4):469–81.
 40
Michalik M, Djahanshiri B, Leo JC, Linke D. Reverse Vaccinology: The Pathway from Genomes and Epitope Predictions to Tailored Recombinant Vaccines. Methods Mol Biol. 2016; 1403:87–106.
 41
Andreoni F, Amagliani G, Magnani M. Selection of vaccine candidates for fish pasteurellosis using reverse vaccinology and an in vitro screening approach. Methods Mol Biol. 2016; 1404:181–92.
 42
Yang YT, Chow YH, Hsiao KN, Hu KC, Chiang JR, Wu SC, Chong P, Liu CC. Development of a fulllength cDNAderived enterovirus A71 vaccine candidate using reverse genetics technology. Antivir Res. 2016; 132:225–32.
 43
Meunier M, GuyardNicodeme M, Hirchaud E, Parra A, Chemaly M, Dory D. Identification of novel vaccine candidates against campylobacter through reverse vaccinology. J Immunol Res. 2016; 2016:5715790.
 44
De Groot AS, Bosma A, Chinai N, Frost J, Jesdale BM, Gonzalez MA, Martin W, SaintAubin C. From genome to vaccine: in silico predictions, ex vivo verification. Vaccine. 2001; 19(31):4385–95.
 45
Lafuente EM, Reche PA. Prediction of MHCpeptide binding: a systematic and comprehensive overview. Curr Pharm Des. 2009; 15(28):3209–20.
 46
GomezBombarelli R, AguileraIparraguirre J, Hirzel TD, Duvenaud D, Maclaurin D, BloodForsythe MA, Chae HS, Einzinger M, Ha DG, Wu T, Markopoulos G, Jeon S, Kang H, Miyazaki H, Numata M, Kim S, Huang W, Hong SI, Baldo M, Adams RP, AspuruGuzik A. Design of efficient molecular organic lightemitting diodes by a highthroughput virtual screening and experimental approach. Nat Mater. 2016; 15(10):1120–1127.
 47
Geanes AR, Cho HP, Nance KD, McGowan KM, Conn PJ, Jones CK, Meiler J, Lindsley CW. Ligandbased virtual screen for the discovery of novel M5 inhibitor chemotypes. Bioorg Med Chem Lett. 2016; 26(18):4487–4491.
 48
Xu Y, Yue L, Wang Y, Xing J, Chen Z, Shi Z, Liu R, Liu YC, Luo X, Jiang H, Chen K, Luo C, Zheng M. Discovery of Novel Inhibitors Targeting MeninMixed Lineage Leukemia (MLL) Interface Using Pharmacophore and DockingBased Virtual Screening. J Chem Inf Model. 2016; 56(9):1847–1855.
 49
Khalili S, Mohammadpour H, Shokrollahi Barough M, Kokhaei P. ILP2 modeling and virtual screening of an FDAapproved library:a possible anticancer therapy. Turk J Med Sci. 2016; 46(4):1135–43.
 50
Bonin CRB, Fernandes GC, dos Santos RW, Lobosco M. A simplified mathematicalcomputational model of the immune response to the yellow fever vaccine. In: Proceedings of the IEEE International Conference on Bioinformatics and Biomedicine. Kansas City: IEEE: 2017. p. 1–8.
 51
Le D, Miller JD, Ganusov VV. Mathematical modeling provides kinetic details of the human immune response to vaccination. Front Cellular Infect Microbiol. 2015; 4:177.
 52
Goutelle S, Maurin M, Rougier F, Barbaut X, Bourguignon L, Ducher M, Maire P. The Hill equation: a review of its capabilities in pharmacological modelling. Fundam Clin Pharmacol. 2008; 22(6):633–48.
 53
Haefner JW. Modeling Biological Systems:Principles and Applications, 1st edn. London: Chapman & Hall, Ltd; 1996.
 54
Bonin C, dos Santos RW, Fernandes G, Lobosco M. Computational modeling of the immune response to yellow fever. J Comput Appl Math. 2016; 295:127–38.
 55
Odeint. Odeint’s Homepage. 2014. http://docs.scipy.org. Accessed 13 May 2018.
 56
LeVeque RJ. Finite Difference Methods for Ordinary and Partial Differential Equations  Steadystate and Timedependent Problems. USA: SIAM; 2007, p. 1341.
 57
Reinhardt B, Jaspert R, Niedrig M, Kostner C, L’ageStehr J. Development of viremia and humoral and cellular parameters of immune activation after vaccination with yellow fever virus strain 17d: a model of human flavivirus infection. J Med Virol. 1998; 56(2):159–67.
 58
Collaborative Group for Studies on Yellow Fever Vaccines. Duration of postvaccination immunity against yellow fever in adults. Vaccine. 2014; 32(39):4977–84.
 59
Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S. Global Sensitivity Analysis: the Primer. New York: Wiley; 2008.
 60
Kay A, Chen LH, Sisti M, Monath TP. Yellow fever vaccine seroconversion in travelers. Am J Trop Med Hyg. 2011; 85(4):748–9.
 61
Akondy RS, Monson ND, Miller JD, Edupuganti S, Teuwen D, Wu H, Quyyumi F, Garg S, Altman JD, Del Rio C, et al. The yellow fever virus vaccine induces a broad and polyfunctional human memory cd8+ t cell response. J Immunol. 2009; 183(12):7919–30.
 62
Monath TP, Gershman M, Staples JE, Barrett ADT. Vaccines (Sixth Edition), 6th. edition In: Plotkin SA, Orenstein WA, Offit PA, editors. London: W.B. Saunders: 2013. p. 870–968.
 63
Monath TP. Yellow fever vaccine. Expert Rev Vaccines. 2005; 4(4):553–74.
 64
OMS. Vaccines and vaccination against yellow fever: Who position paper–june 2013. Wkly Epidemiol Rec. 2013; 88:269–84.
 65
de Souza Lopes O, de Almeida Guimarães SSD, de Carvalho R. Studies on yellow fever vaccine iii—dose response in volunteers. J Biol Stand. 1988; 16(2):77–82.
 66
Vieira P, Rajewsky K. The bulk of endogenously produced igg2a is eliminated from the serum of adult c57bl/6 mice with a halflife of 6–8 days. Eur J Immunol. 1986; 16(7):871–4.
 67
Vieira P, Rajewsky K. The halflives of serum immunoglobulins in adult mice. Eur J Immunol. 1988; 18(2):313–6.
 68
Lee HY, Topham DJ, Park SY, Hollenbaugh J, Treanor J, Mosmann TR, Jin X, Ward BM, Miao H, HoldenWiltse J, Perelson AS, Zand M, Wu H. Simulation and prediction of the adaptive immune response to influenza A virus infection. J Virol. 2009; 83(14):7151–65.
Funding
Authors would like to thank the support from CNPq, FAPEMIG, and UFJF.
Availability of data and materials
Data sharing not applicable to this article as no datasets were generated or analysed during the current study. Model parameters and initial conditions used in simulations are included in this published article. Code to solve the mathematical model can be made available upon request to the authors.
Author information
Affiliations
Contributions
Conception and design of the mathematical model: all authors. Computational implementation of the mathematical model: CRBB. CRBB has also been involved in drafting the manuscript. All authors read, revised and approved the final manuscript.
Corresponding author
Correspondence to Carla R. B. Bonin.
Ethics declarations
Ethics approval and consent to participate
Not applicable
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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.
About this article
Received
Accepted
Published
DOI
Keywords
 Computational vaccinology
 Yellow fever
 Mathematical modeling
 Computational modeling
 Immune system
 Ordinary differential equations