Understanding Original Antigenic Sin in Influenza with a Dynamical System

Original antigenic sin is the phenomenon in which prior exposure to an antigen leads to a subsequent suboptimal immune response to a related antigen. Immune memory normally allows for an improved and rapid response to antigens previously seen and is the mechanism by which vaccination works. I here develop a dynamical system model of the mechanism of original antigenic sin in influenza, clarifying and explaining the detailed spin-glass treatment of original antigenic sin. The dynamical system describes the viral load, the quantities of healthy and infected epithelial cells, the concentrations of naïve and memory antibodies, and the affinities of naïve and memory antibodies. I give explicit correspondences between the microscopic variables of the spin-glass model and those of the present dynamical system model. The dynamical system model reproduces the phenomenon of original antigenic sin and describes how a competition between different types of B cells compromises the overall effect of immune response. I illustrate the competition between the naïve and the memory antibodies as a function of the antigenic distance between the initial and subsequent antigens. The suboptimal immune response caused by original antigenic sin is observed when the host is exposed to an antigen which has intermediate antigenic distance to a second antigen previously recognized by the host's immune system.


Introduction
An immune memory comes from a previous infection or vaccination, stores the information for antigen recognition, and is activated in a future infection by a related pathogen. Long-term immune memory has been observed in various pathogens including smallpox [1], malaria [2], hepatitis B [3], dengue [4], and influenza A [5]. By recognizing and rapidly eliminating the reinfecting pathogen, this long-lasting effect can permanently or temporarily prevent the reinfection of the host by some pathogens [6]. In some cases, this long-lasting effect can also reduce the severity, duration, or risk of the infection and symptoms [7]. Smallpox virus, also called variola virus, only propagates in humans and has a relatively low mutation rate [8]. In contrast, influenza A virus propagates in humans, pigs, and aquatic birds, with a higher mutation rate that is approximately 2:0|10 {6 / nucleotide/infectious cycle [9], or 1:6|10 {5 /amino acid/day. Calculation of the binding free energy between human antibodies and circulating influenza A strains shows that the virus mutates away from the genotypes that code for hemagglutinin proteins well recognized by the human immune system [10]. Thus for influenza A, there is usually a significant antigenic distance between the circulating strain in a given year and the immune memory from previous years.
Original antigenic sin is the phenomenon in which prior exposure to an antigen leads to a subsequent suboptimal immune response to a related antigen [11][12][13]. In some years when the antigenic distances between vaccine and circulating virus strains fell into a certain range, the effect of original antigenic sin decreased the effectiveness of influenza vaccines. Historical data of influenza vaccines indicate that vaccine effectiveness does not monotonically decrease with the antigenic distance between the vaccine strains and the circulating strains, but rather has a minimum at an intermediate antigenic distance [14,15]. Interestingly, since the vaccine effectiveness at this intermediate antigenic distance between the vaccine and circulating strains is lower than the effectiveness at a larger antigenic distance in unvaccinated people, original antigenic sin could make vaccinated people more susceptible to the virus than those who are unvaccinated.
The mechanism of original antigenic sin was previously studied using stochastic models at the cellular level [16,17]. These previous studies developed stochastic models with thousands to millions of B cells [16,17]. The stochastic models introduce various antigens to a repertoire of B cells. The B cells with higher affinity to an antigen have larger probability to be selected during the B cell maturation process. Earlier works discussing the mechanism of original antigenic sin at the cellular level include [17], which attributed original antigenic sin to the localization of the B cells in the secondary immune response around the B cells in the primary immune response in the amino acid sequence space. The affinity between an antibody and an antigen is given by the generalized NK model (GNK model) of the three-dimensional protein structures [18]. The GNK model was derived from the NK model which was originally introduced to model rugged fitness landscapes [19,20] and evolutionary processes [21][22][23]. In the GNK model, the amino acid sequences of a group of influenza A specific antibodies are allowed to mutate freely and independently in the affinity landscape to maximize their affinities to the virus. B cells that produce antibodies with the highest affinities replicate into the next generation. The mutation of the virus is modeled by changing the fitness landscape. The antibody affinities at the end of the simulation correlate well with the vaccine effectiveness data observed in history [14].
The present study aims to use a set of ordinary differential equations (ODEs) to describe the interaction among the B cells, the virus particles, and the epithelial cells. This deterministic model can reduce the memory and CPU requirements, compared to the stochastic models [16,17]. In this study, I developed a set of ODEs as the mean-field approximation of the stochastic models that store the information of each B cell. The ODEs give a deterministic explanation of original antigenic sin. This explanation agrees with both the observed data [14] and a previous explanation given by the GNK model [17]. Various ODE models have been established to describe and simulate the process of influenza A infection and the resulting immune response [24][25][26][27][28][29]. The basic elements of these ODE models were described in [30]. Automata have also been used to model the time-dependent spatial distribution of the tissue cells and the virus [31,32].
In this study, I build a deterministic ODE-based model compatible with the previous GNK model [17] to reproduce the observed phenomenon of original antigenic sin. The main purpose of this paper is to address the following two questions on the modeling of original antigenic sin. First, can a deterministic dynamical model with a small number of state variables simulate the phenomenon of original antigenic sin, which was simulated by stochastic models of a large repertoire of B cells [16,17]? Second, what is the mechanism of original antigenic sin revealed by the deterministic dynamical model? I address the first question by building a deterministic ODE model, which reproduces the phenomenon of original antigenic sin observed in experiment [11][12][13]. See the subsections Model Development and Description, Time Courses of Infection and Recovery, and A General Picture of Original Antigenic Sin in the section Results. I address the second question by analyzing the non-monotonicity of the overall effect of an immune response. See the subsection Mechanism of Original Antigenic Sin in the section Results. The values of the parameters mainly come from previous studies that simulated influenza infection and immune response and obtained plausible results. However, the limitations of the available experimental data do not allow one to develop an accurate model purely based upon experimental data [33]. Experimental data available for model development usually have limited quality or quantity, causing unavoidable overparameterization of the models. In the present dynamical model, two parameters c 0 and s need to be estimated prior to the simulation. I give the estimation of parameters c 0 and s in the section Materials and Methods and perform the sensitivity analysis of both of them in the section Results. The terms in the ODEs have clear physical meanings, so my model explicitly illustrates the details of the influenza A infection and the immune response. A comparison between this deterministic dynamical model and previous stochastic models is presented in the section Discussion. A brief review of the influenza A genome is in Appendix S1.

Model Development and Description
I use a simplified model consisting of the major components of an immune response, which are epithelial cells, influenza A viruses, and an immune system, to describe the dynamics of an influenza A infection and the subsequent immune response. This model contains six state variables, which are the healthy cell concentration H ð Þ, the infected cell concentration I ð Þ, the viral load V ð Þ, the concentrations of naïve and memory antibodies X 1 and X 2 ð Þ , respectively, and the naïve antibody affinity U 1 ð Þ. These state variables are in Table 1.
This model comes from the following information on influenza A infection. The concentration of epithelial cells on the upper respiratory tract is around a fixed homeostatic level H 0 , which is also the sum of the concentrations of healthy cells (H), of infected cells (I), and of dead cells (D) killed by the influenza A virus, respectively. Free influenza A virus particles (V ) are released from infected cells and are eliminated by influenza A specific antibodies. I only consider virus clearance by antibodies, because a T cell can recognize influenza A strains with different antigenic characters [34]. This model has two types of antibodies: the naïve antibodies and the memory antibodies generated by the last influenza A infection or vaccination. The concentrations of the naïve and the memory antibodies are defined as X 1 and X 2 , respectively. The immune system cleared the influenza A viruses bound by antibodies. With the definition of antibody affinity the concentration of influenza A virus particles bound by antibodies Ag : Ab ½ is proportional to the concentrations of the free influenza A virus particles Ag ½ and of the influenza A specific antibodies Ab ½ . The naïve and memory antibodies have affinities U 1 and U 2 , respectively, to the influenza A virus. The affinity of memory antibodies, U 2 , is a constant parameter of the model. The maximum affinity is defined as U max . Here U 1 and U 2 are quantified using the reduced unit 1:0|10 7 M {1 , as described in the subsection Reduced Units and Parameter Estimation in the section Materials and Methods. Thus the affinities U 1 and U 2 are defined as U i~Ka =1:0|10 7 M {1 (i~1,2). The maximum affinity is U max~1 when K a~1 :0|10 7 M {1 .
From the above information, a minimal set of ODEs are built to model the influenza A specific immune response of co-existing naïve antibodies with a low initial affinity and memory antibodies with a higher and constant affinity. The state variables Z~H,I,V ,X 1 ,X 2 ,U 1 ð Þcomprise the healthy cell concentration H, the infected cell concentration I, the viral load V , the naïve antibody concentration X 1 , the memory antibody concentration X 2 , and the binding affinity of naïve antibodies U 1 .
The homeostasis of epithelial cells gives an additional algebraic equation for the dead cell concentration, D, Equation 2 describes the dynamics of the concentration of healthy epithelial cells. The repair mechanism for epithelial tissue is activated only if any damage in epithelial cells is detected (D=0), and new healthy cells are regenerated with the rate lD to repair the damaged tissue [25]. Alternative models for this repair mechanism include regeneration rate~l [28] or lDH [29]. In the stochastic model in [32], the regeneration rate is 0 when D~0 and has a mathematical expectation of lH when D=0. The average life span of human trachea cells is 47.5 days [35], while the time for epithelial cell regeneration is 0.3-1 day [25,36], showing that the cell regeneration rate is significantly higher than the normal cell division rate. With this consideration in mind, I select the expression lD as the cell regeneration rate [25]. The infection rate b represents the speed in which influenza A virus converts healthy cells into infected cell. The protective effect of interferon is neglected in this simplified model. Equation 3 characterizes the dynamics of the infected cell concentration. All the infected cells are converted from healthy cells. Infected cells are killed by the virus with the rate a.
Equation 4 depicts the generation and elimination of virus particles. Virus particles are released from infected cells with the rate k. The half-life of free virus particles is 1=m. Viruses bound by antibodies are neutralized and cleared by the immune system. Thus the virus clearance rate is proportional to the concentration of viruses bound by antibodies Ag : Ab ½ . From equation 1, the virus clearance rate is proportional to the antibody affinity U i (i~1,2), the antibody concentration X i (i~1,2), and the viral load V , respectively.
Equations 5 and 6 show the secretion and decay of naïve antibodies and memory antibodies. Antigen presenting cells (APC) process the virus and present the antigen on their surface, activating naïve T cells. Some of these activated T cells proliferate and differentiate into Th2 helper T cells. Th2 cells and free virions activate B cells together [37]. The intensity of activation signal for B cells, c Z,t ð Þ, is a function of time depending on viruses, APCs, and naïve T cells. The intensity c Z,t ð Þ is a rectangular window function with the maximum value of c 0 , and is further described in the section Materials and Methods. Naïve B cells mature in germinal centers, undergoing proliferation and somatic hypermutation. B cells are selected by competing for antigen binding and activation signals from Th2 cells surrounding the germinal center. The morphology of germinal centers determines that the interface between the B cell region and the Th2 cell region is approximately constant, and so is the amount of antigens inside the germinal center. Therefore B cells inside the germinal center compete with each other for the activation signal. The ratio of the intensities of the activation signals for naïve and for memory B cells is The decay rate of both naïve and memory antibodies is b. I use identical decay rate (b) for naïve antibodies (X 1 ) in equation 5 and memory antibodies (X 2 ) in equation 6 because the decay rate is independent of the type of antibodies [38]. Equation 7 indicates that the increase rate of the affinity is proportional to the concentration of the antigen-antibody complex. Because the B cells are selected by the affinity to the antigen in their maturation process, the increase of the naïve antibody affinity U 1 is driven by successful binding between the naïve antibody and the antigen. The logistic factor U max {U 1 ð Þ ensures that the probability for B cells to mutate to a state of higher affinity decreases as the maturation proceeds.
The dynamical model comprises equations 2-7. Equations 2, 3, 4, and 7 are adapted from previous models [24,25,29]. Note that equations 2, 3, 4, and 7 are not identical to their original form in literature. Most of these previous models were developed according to the processes of cell death and regeneration and virus entry and release. The parameters of these models come from experiment except for the parameter c 0 , which describes the activation of B cells, and for the parameter s, which describes the B cell maturation process. The parameter s was fit to experimental data in a previous model [29]. As will be shown, I estimate the values of c 0 and s all over again and perform a sensitivity analysis to the parameters c 0 and s.
The present six-ODE model is able to be mapped from a previous 10-ODE model developed by Hancioglu et al.

Time Courses of Infection and Recovery
With all parameters defined and fixed in the section Materials and Methods, I use the stiff differential equation solver ode23s in MATLAB to numerically solve equations 2-7. The relative and absolute error tolerances of the solver are 10 {9 and 10 {6 , respectively. The first set of parameters listed in the left column of Table 2 I run a simulation of 20 days. The solved trajectories of the state variables Z are compared to the kinetics of influenza A infection observed in reality to verify the model parameters. I use two cases to illustrate the dynamics of all state variables. In the first case, the virus has substantially mutated from the previous strains, and the binding affinity of the memory antibodies to the virus is low. In the second case, there is no significant escape mutation of the virus, and the binding affinity of the memory antibodies to the virus is high. The affinity of memory antibodies is U 2~1 0 {3 in the first case and is U 2~0 :5 in the second case. The details of the model dynamics are shown in Figure 1 and 2. Figure 1 describes the whole process of influenza A virus infection and clearance in humans without immune memory (U 2~1 0 {3 ). A symptom with approximately 30% of the epithelial cells killed is observed after the infection. The peak of the dead cell concentration D occurs on Day 1 to Day 2, agreeing with the experimental data [25]. On Day 5, the dead cell concentration D falls under 0.1. The viral load V decreases to the initial level on Day 3 to Day 4. A 10 5 -fold increase of the naïve antibody concentration occurs in the virus infection and clearance process. The naïve antibody affinity U 1 approaches to the maximum U max~1 . The memory antibody concentration X 2 has an initial 10 2 -fold increase, and decreases approximately exponentially after Day 1 with the rate 0.0427, similar to the antibody decay rate b~0:043. Thus few memory antibodies are produced after Day 1. Figure 2 shows the dynamics of virus infection and clearance with immune memory (U 2~0 :5). Few dead cells are accumulated and thus no symptoms are observable in the infected person. The viral load is remarkably suppressed compared to Figure 1. Figure 2 depicts the effect of a successful vaccination. Compared to Figure 1, the increase of the naïve antibodies concentration is absent, and the naïve antibody concentration, X 1 , decreases approximately exponentially with the rate 0.0423, close to the antibody decay rate b~0:043, indicating that naïve antibodies are rarely produced during the whole process of virus infection and clearance. No significant somatic hypermutation is observed in those naïve B cells. The naïve antibody affinity U 1 is almost constant, as shown by the plot of U 1 against t. There is a notable increase in the memory antibody concentration X 2 : the value of X 2 on Day 7 is approximately 10-fold higher than that in Figure 1.
The immune response is dominated by the naïve antibodies when the memory antibodies have low affinity and is dominated by the memory antibodies when they have high affinity. The transition between these two cases occurs when the value of the memory antibody affinity U 2 falls into a critical region. In the following subsection, U 2 is set to a variety of values in the range 10 {3 ƒU 2 ƒ1. As the value of U 2 changes, the phenomenon of original antigenic sin can be observed in other characters of the dynamics, such as the maximum percentage of dead cells, the maximum viral load, the cumulative effects of naïve antibodies and of memory antibodies, and the average antibody affinity. This model is able to reproduce the phenomenon of original antigenic sin observed in the experimental data at intermediate memory antibody affinity U 2 .

A General Picture of Original Antigenic Sin
To illustrate the phenomenon of original antigenic sin, I choose 100 values of U 2 logarithmically spaced between 10 {3 and 1:0. The minimum value U min 2~1 0 {3 reflects the case that memory antibodies rarely recognize a new virus strain. The maximum value U max 2~1 :0 corresponds to the immune response with the highest memory antibody affinity. The intermediate values of U 2 correspond to the case that memory antibodies have decreased capability to recognize a new virus strain due to the escape mutation of the strain or imperfect vaccination. One hundred independent simulations were run with these 100 values of U 2 , respectively. The maximum viral load and the maximum percentage of dead cells were recorded for each simulation. The cumulative effects of naïve antibodies and of memory antibodies are respectively calculated with Similarly, the average antibody affinity at the end of each simulation is Equation 11 calculates the average antibody affinity after the 20day period of time during which the patient recovered from the infection. Note that indicated by equation 7, the naïve antibody affinity U 1 monotonically increases, while the memory antibody affinity U 2 is constant. Figure 3 depicts the maximum percentage of dead cells, the maximum viral load, the cumulative effects of naïve and of memory antibodies, and the final average antibody affinity as the functions of the memory antibody affinity U 2 . The peak in the figure, 40.8%, is reached when U 2~3 :8|10 {3 and is 33.3% higher than the maximum percentage D max 0~3 0:6% with U 2~1 0 {3 . This 33.3% increase is significant for the percentage of dead epithelial cells. When U 2 w0:1, the maximum percentage of dead cells falls below 10%, indicating that no observable symptoms occur in the infected person. That is, when the antigenic distance from the previous infection or vaccination to the new virus strain is small, the immune system can clear the virus effectively. Figure 3b shows the non-monotonicity of the maximum viral load V as a function of the memory antibody affinity U 2 in the process of virus infection and clearance. The maximum viral load is V max 0~4 5:0 when U min 2~1 0 {3 . The maximum viral load in the region 10 {3 vU 2 v1:2|10 {2 is higher than V max 0~4 5:0, and so original antigenic sin occurs in this region. The maximum viral load in the interval 1:6|10 {3 vU 2 v2:8|10 {3 is at least twice as high as V max 0 . With U 2 w0:11, the maximum viral load is less than unity, agreeing with the fact that memory antibodies effectively recognize and eliminate virus strains which are antigenically similar to the immune memory [14]. Figure 3c describes the cumulative effects, X Int 1 and X Int 2 , of naïve and memory antibodies in the process of virus infection and clearance, respectively. The cumulative effects X Int 1 and X Int 2 are calculated with equations 9 and 10, respectively. There is a   Figure 3c, this sudden decrease in the cumulative effect of naïve antibodies X Int 1 occurs in the same region of U 2 . When U 2 w8:7|10 {3 , U avg increases approximately linearly with U 2 , which is mainly due to the contribution of the memory antibodies.

Mechanism of Original Antigenic Sin
The dynamical system defined by equations 2-7 can be split into two subsystems, i.e. an actuator and a controller, with weak coupling between them. Equations 2-4 constitute the actuator with H, I, and V as the state variables. Equations 5-7 constitute the controller with X 1 , X 2 , and U 1 as the state variables. The actuator is controlled by the variable E~U 1 X 1 zU 2 X 2 w0. Therefore equation 4 is equivalent to The actuator consisting of equations 2, 3, and 12 has two steady states: By calculating the eigenvalues of the Jacobian of the actuator, I find that the first steady state is stable for any Ew0, and the  Figure 4 displays the effect of E: all the viruses are cleared from the patient when E is larger than 0.18, and the decay rates of the dead cell concentration D and of the viral load V increase with E. Therefore, E is the only factor controlling the values of D and V .
The other subsystem is the controller comprising the state variables X 1 , X 2 , and U 1 . The controller observes the state of the actuator as the factor c Z,t ð Þ, which jumps from 0 to c 0 when the viral load V reaches 0.1 and remains c 0 for 14 days. Due to the quick virus proliferation at the beginning of the infection, I let c Z,t ð Þ~c 0~1 as an approximation. The dynamics of the expression E are described by the following equation: The first term on the right hand side of equation 13, U 2 d X 1 zX 2 ð Þ =dt, is the product of U 2 and the derivative of a first order process X 1 zX 2 . The derivative of X 1 zX 2 is obtained by adding equation 5 to equation 6 with the approximation c Z,t ð Þ~1 and is independent of U 2 . The form of equations 5 and 7 shows that U 2 suppresses the variables X 1 t ð Þ and U 1 t ð Þ, and so the term X 1 dU 1 =dt monotonically decreases with U 2 . In the case of small U 2 , the factor U 1 {U 2 ð Þ w0 in the process of virus infection and clearance (see Figure 1), thus the third term U 1 {U 2 ð Þ dX 1 =dt decreases with U 2 . In the case of large U 2 , X 1 is approximately constant, and the third term is negligible. Consequently, the first term in equation 13 increases with U 2 and the other terms decrease with U 2 . Figure 5 shows E as a The source of original antigenic sin is the interaction of the state variables of the controller, or the immune system, modeled by equations 5-7. When the memory antibody affinity U 2 is either low or high, naïve or memory antibodies are dominant, respectively. If naïve antibodies are dominant, their final affinity is high. At intermediate U 2 , the interaction and competition between naïve and memory antibodies lead to a decreased immune effect E, which clears the viruses less effectively. This is original antigenic sin. When original antigenic sin occurs, the influenza illness rate increases [14] due to the increases of D and V during the process, and the average antibody affinity decreases [17] due to the decrease of E. Note that the average antibody affinity is defined as U avg~E = X 1 zX 2 ð Þ , where X 1 zX 2 ð Þ are approximately independent of U 2 , as discussed above.

Sensitivity Analysis
The parameters other than c 0 and s come from literature. See Table 2. As mentioned in the subsection Model Development and Description, my model described by equations 2-7 can be mapped from a previous dynamical model of influenza infection and immune response [29], in which a comprehensive sensitivity analysis for parameters is available. The sensitivity analysis for most of the parameters in this model has been performed in [29]. So I first focus on two remaining parameters, c 0 and s. The parameter c 0 characterizes the stimulation of the immune system when the viral load increases beyond a threshold. Because this paper aims to give a simplified model comprising the most important factors of both epithelial cells and an immune system, the effects of APC and Th2 cells are combined into the parameter c 0 . The parameter s reflects the process of B cell somatic hypermutation which produces antibodies with high affinity to the antigens. Here I present a sensitivity analysis for the parameters. Figure 6 describes the behavior of the dynamical system with different values of parameters c 0 and s. The maximum percentage of dead cells with large U 2 is insensitive to s. The average antibody affinity with small and large U 2 is also insensitive to both c 0 and s, but is sensitive to intermediate U 2 . If U 2 is higher than a threshold, the memory antibodies play the major role in the immune response; otherwise the naïve antibodies mainly conduct the immune response. As shown in Figure 6b and 6d, this threshold of U 2 decreases with c 0 and increases with s. However, the dynamics with different values of c 0 and s in Figure 6 resemble those in Figure 3a and 3d. The existence of original antigenic sin is insensitive to the parameters c 0 and s. In addition, Figure 7 shows that both the maximum percentage of dead cells and the average antibody affinity U avg are insensitive to the decay rate of antibodies (b) and the initial concentration of memory antibodies (X 2 0 ð Þ). This sensitivity analysis shows that the severity of an influenza A infection decreases with c 0 , and the effect of original antigenic sin increases with c 0 and decreases with s.

Discussion
The model defined by equations 2-7 is a significant simplification to the previous models describing the kinetics of influenza A infection [25,29]. This model introduces a second type of antibodies, the memory antibodies, to simulate the competition and cooperation between naïve and memory antibodies. The contributions of Th1 cells, cytotoxic T lymphocytes (CTLs), interferons, and epithelial cells protected by interferons are incorporated into the parameters of the model. The contributions of APCs and Th2 cells, which together activate B cells, are captured by the factor c Z,t ð Þ in a mean-field approach, rather than being explicitly modeled using ODEs. The presenting model has two limitations. First, at the same time of the increase in the dead cell concentration D, the viral load V increases by 10 3 -10 4 fold and reaches the maximum on Day 1, which is different from the experimental results of 10 4 -10 5 fold increase at the maximum on Day 2. Second, this model does not contain a term modeling the loss of antibodies due to virus binding in equations 5 and 6. These limitations are due to simplification of the model. However, these limitations do not seriously affect the emergence of original antigenic sin in the model at intermediate memory antibody affinity U 2 , the major topic of this study.
The concentration of CTLs increases by 100 times in the first seven days after infection [25] to eliminate the infected cells. The cellular immune system usually has strong cross immunity for antigenically different influenza A virus strains [25], while the humoral immune system cannot effectively recognize a new influenza A virus strain with a large antigenic distance from the previous strains seen by the host immune system [14]. Thus the effects of CTLs against different influenza A virus strains are more homogeneous than those of antibodies. That is, CTLs induced by previous virus strains can effectively suppress a new virus strains despite the escape mutation of the virus [25]. By contrast, the antibody affinity decreases with the antigenic distance between the previous virus strains and the new strain. Thus the contribution of CTLs is more constant compared to that of antibody and can be modeled by constant parameters to describe original antigenic sin.
For the same reason, the protection for healthy cells by the interferon secreted by infected epithelial cells is not modeled as an independent equation either. Additionally, the interferon and the cells protected by the interferon are not the key factors of the dynamics of virus infection and clearance: an absence of interferons does not affect the final elimination of all viruses and dead cells [29]. APCs and Th2 cells have little interaction with the elements in the model other than the antibodies. Hence a simple function c Z,t ð Þ is introduced to model the activation of B cells. The present model contains two parts. Equations 2-4 constitute a general model for an infection in tissue caused by a cytopathic virus. Equations 5-7 define a model for the immune system which recognizes and clears the virus. This model can be extended to include the dynamics of CTLs and interferons in the immune system. The model can also be extended to simultaneously consider multiple types of cytopathic viruses by using equations 2-4 with different parameter sets l,b,a,k,m ð Þto model each type of virus. The cytopathic viruses fall into two categories: those causing acute diseases and those causing chronic diseases. Viruses in the first category, such as influenza A virus, are cleared in a short period of time, thus relatively few escape mutations occur. The immune system therefore directs itself towards a fixed target. Equations 5-7 modeling the immune system do not require any modification to take into account the escape mutation of the virus. On the other hand, viruses in the second category, including HIV, persist for years in the host and keep mutating away from the immune system. A new set of ODEs are needed to model this case. First, additional terms are required for equation 7 to describe the decrease of the memory antibody affinity due to the escape mutation of the virus. Second, if the immune system is also infected by the virus, equations 5 and 6 should also contain terms to model this infection. A similar model of the immune response against HIV and the competition between antibodies has been developed [39].
The currently available mathematical models of original antigenic sin falls into two categories: stochastic models represented by the GNK model [17], and deterministic models as the present one. Now I compare the mathematical form of the GNK model [17] with that of my deterministic model. Both models consider the contributions of naïve and memory B cells. Both models explicitly simulate the competition between different types of B cells. In the deterministic model, the maturation of naïve B cells follows a logistic process. In the GNK model, the B cells have random walks on a rugged and random landscape [17,19] where the density of neighboring states with higher affinities decreases with the affinity in the current state. The deterministic model has two variables, U 1 and U 2 , for the naïve and memory antibody affinities, respectively. The GNK model, however, stores the amino acid sequences of 1000 naïve antibodies and other 1000 memory antibodies [17]. After a simulation of 30 generations, the number of different amino acid sequences generally converges to less than five, similar to the deterministic model in which both naïve and memory antibodies are considered as monoclonal.
There are parallels and differences between my deterministic model and the GNK model [17], although these two models have different mathematical forms as shown above. First, the deterministic model uses the factors U 1 X 1 = U 1 X 1 zU 2 X 2 ð Þ and U 2 X 2 = U 1 X 1 zU 2 X 2 ð Þto model the selection of naïve and memory B cells, respectively. In an immune response, B cells producing antibodies with high binding affinity U 1 or U 2 are selected. As a comparison, the stochastic model stores the amino acid sequence of each B cell, selects those B cells with high affinity to the antigen, and replicates the selected B cells to the next generation. Second, the deterministic model simulates the B cell maturation process with the factor c Z,t ð Þ describing the activation of B cells. In the simulation, c Z,t ð Þis greater than zero for 14 days, in which the B cells compete for the antigen and divide. The stochastic model repeats the process of B cell hypermutation and selection for 30 generations of B cells, which correspond to the primary or secondary immune response. Third, the naïve antibody affinity in the deterministic model, U 1 , is modeled by equation 7, a logistic equation. The increase rate of U 1 decreases as U 1 approaches to the maximum antibody affinity, U max . The stochastic model builds a rugged antibody affinity landscape, in which the locations with high affinities have low density. Consequently, the deterministic model is a mean-field approximation of the B cell maturation process. The deterministic model is able to simulate original antigenic sin with reduced memory and CPU requirement, while ignoring the amino acid sequence of each B cell.
The dynamical model introduced in this paper is deterministic, while the process of influenza A infection and clearance is stochastic in nature. However, the deterministic model gives similar results as the stochastic models [16,17]. The deterministic model assumes both naïve and memory antibodies to be monoclonal. If either naïve or memory antibodies have multiple amino acid sequences, the ODEs could be modified by introducing more types of antibodies. Hence the competition factors U 1 X 1 =(U 1 X 1 zU 2 X 2 ) and U 2 X 2 =(U 1 X 1 zU 2 X 2 ) could be replaced by a tournament-like algorithm involving all the antibodies recognizing the antigen. By introducing the method of splitting the dynamical system into an actuator and a controller, the present model provides a starting point for the application of nonlinear control theory to explain original antigenic sin. The dynamical model can also be helpful to rational vaccine design.

Kinetics of Influenza A Virus Infection
Influenza A virus infection can be described by a dynamical process. The infection occurs in the epithelial cells on the surface of upper respiratory tract in the bronchi with diameter larger than 3.3 mm [25]. The incubation period between the infection and the emergence of symptoms ranges from one day to five days and is typically two days. The host starts to shed infectious virus particle approximately 24 hours prior to the emergence of symptoms. The typical initial concentration of influenza A virus particles is 10 {13 M. The viral load usually reaches a maximum 3|10 {9 M two days after infection and falls back to the initial level six days after infection [25]. Influenza A virus is cytopathic and kills the infected cells, causing the dead cells to accumulate in situ. The percentage of dead epithelial cells reaches the maximum of 30-50% on Day 2 and decreases to 10% on Day 5. If the maximum percentage is lower than 10% in the process of virus infection and clearance, no symptoms are observable [25]. The immune system is activated by the detection of virus particles. The most important suppressor of influenza A virus are antibodies IgG and IgA, followed by the CD8 CTLs [37]. The concentrations of B cells and of plasma cells increase by 10 2 times and 2|10 4 times, respectively, within seven days [25]. The immune response to a primary infection generates memory antibodies with binding affinity 10 6 M {1 and concentration 10 {13 M, constituting 0.1%-1% of the total antibodies [29,37]. The naïve antibodies capable of recognizing the antigen have the affinity 10 4 M {1 and concentration 10 {15 M [29,37], constituting 0.001%-0.01% of the total antibodies.

Reduced Units and Parameter Estimation
I use reduced units for all the variables and parameters to make their values close to unity and to facilitate the numerical calculation. For the state variables H, I, V , X 1 , and X 2 , the unit is defined as the homeostatic concentration of epithelial cells in the upper respiratory tract, which is 1:7|10 {11 M [25,29]. The unit of U 1 and U 2 is defined as the maximum affinity between memory antibodies and influenza A viruses, which is 1:0|10 7 M {1 [37]. The reduced units for all the variables in equations 2-7 are listed in Table 1.
The majority of the parameters in this dynamical model are obtained from previous experiments. These parameters fall into two sets that are compatible with each other. The first set of parameters were given by the publications [24,25,29,36]. These publications depicted the process of influenza A virus infection and clearance in the cellular level, taking into account the concentrations of epithelial cells, viruses, APCs, interferons, Th1 and Th2 helper cells, CTLs, B cells, plasma cells, and antibodies. The second set of parameters were extracted from an influenza A virus infection experiment with six volunteers [27]. A simpler ODE model with a fixed parametric form was built to fit the daily viral loads data measured from nasal wash [27]. These two sets of parameters in the reduced units are listed in Table 2. Despite the different approaches to obtain the parameters, the parameters b, a, k, and m from [24,25,29,36] and [27] are similar.
Compared to some previous models [25,29], a major simplification in this study is neglecting the propagation of the activation signal for the immune system, originated by the detection of the virus and through APCs, Th2 cells, and B cells. Instead, I introduce a time-dependent factor c Z,t ð Þ to model the activation signal for the immune system in a mean-field approach. In a typical process of influenza A virus infection, the viral load and the concentration of APCs reach the maximum simultaneously on Day 2 [25]. The viral load decreases to the initial level on Day 6 [25]. As listed in Table 3, the half-lives of APCs, helper T cells, B cells, and plasma cells are similar to the duration of the infection. Thus the duration of B cell maturation process is estimated to be 14 days, longer than the duration of viral clearance. Accordingly, the factor c Z,t ð Þ has the initial value of zero, is assigned the value c 0 when V reaches 0.1, and equals to c 0 for 14 days before decreasing to zero again. Using the output of the previous models [25,29], I estimate the parameter c 0 to be 1.0, and the parameter s to be 100.