Research Article | Open Access

# Immune Response to a Variable Pathogen: A Stochastic Model with Two Interlocked Darwinian Entities

**Academic Editor:**Zvia Agur

#### Abstract

This paper presents the modeling of a host immune system, more precisely the immune effector cell and immune memory cell population, and its interaction with an invading pathogen population. It will tackle two issues of interest; on the one hand, in defining a stochastic model accounting for the inherent nature of organisms in population dynamics, namely multiplication with mutation and selection; on the other hand, in providing a description of pathogens that may vary their antigens through mutations during infection of the host. Unlike most of the literature, which models the dynamics with first-order differential equations, this paper proposes a Galton-Watson type branching process to describe stochastically by whole distributions the population dynamics of pathogens and immune cells. In the first model case, the pathogen of a given type is either eradicated or shows oscillatory chronic response. In the second model case, the pathogen shows variational behavior changing its antigen resulting in a prolonged immune reaction.

#### 1. Introducing Relevant Prior Knowledge

##### 1.1. Putting the Objectives of the Paper into Context

The wide relevance of pathogens, such as the influenza virus, the human immunodeficiency virus (HIV), or trypanosomes, give great significance to those studies, where pathogens are able to vary their antigens while still vital in the host and where the host’s immune system mounts specific immune reactions (by clonal selection, somatic hyper-mutation, and forming an immune memory) [1].

Investigations of long-term dynamics of hosts and their immune systems in environments that consist of variable pathogen strains are especially valuable in, first, knowing how duration of the immunological memory can influence the pathogen competition and in, second, evaluating whether the pathogen can be a selective force that can shape the evolution of the immunological memory [2]. The study of these processes is, however, a very complex endeavor. Indeed, in the lowest approximation of understanding the interaction between the invading pathogen and the immune system, the selected immune clones do not go on to future generations of the infected host. Moreover, the ability of a virus/bacteria to survive within the host does not necessarily imply good ability to infect other hosts, and thus survive and evolve.

In this paper, we will focus solely on modeling the dynamics of an infection within one host, and we will provide possible understanding of how the pathogen load and pathogen diversity influence the immune response [3–5]. Can the complex process of an immune response be simplified to be tractable theoretically but still represent some basic facts from immunobiology [6]? In understanding the immune response, it is well established that both the pathogen [7] invading the host as well as the effector [8, 9] of the host’s immune system (trying to get rid of the pathogen) undergo a step-by-step Darwinian process, namely, multiplication with mutation, and selection. This process is stochastic in nature: chance events weighted by fitness influence the processes of multiplication, mutation and selection. The immune response involves two such entities, which are coupled: the pathogen, that is, virus, bacterium or parasite, on the one hand, and the immune effector cell together with its immune memory cell as idioblast on the other hand. The specific immune response to the pathogen worsens the conditions for the pathogen to thrive, and ultimately eliminates the pathogen, at best, without harming the host.

In the following section (Section 1.2), we provide a short description of the basic biological facts of an immune response as well as some mathematical background on continuous models studied previously in theoretical immunology (Section 1.3). We then propose on grounds of a simple stochastic approach of a Darwinian entity (Sections 2.1–2.4), a stochastic model of an immune response (Section 3.1) by coupling two Darwinian entities. We apply this model to a nonvarying pathogen (Section 3.2), and to the challenging problem of a variable pathogen (Section 3.3), for example, a strain of a pathogen transforming into another strain each with different antigens that are presented to the immune system. Finally we model the maturation process from a naive immune cell to an effector cell that contributes to the elimination of the pathogen (Section 4).

##### 1.2. Basic Facts from Immunology and the Request for a Simple Model

The interaction between a pathogen, which can be a virus, a bacterium or a parasite that has invaded a host, and the reaction of the host’s immune system, which is a concerted action of multiple players in time and space, is certainly not simple [1]. It includes the fully developed specific adaptive/acquired immune system, mainly the B and T lymphocytes as well as the innate immune system, mainly the macrophages, which are dumping cells, and the soluble cytokines, which themselves have a wide spectrum of biological activities that help to coordinate the complex immune regulation.

An important part of the specific adaptive/acquired immune system is the “endogenous-cellular” path, where the pathogen—which is usually a virus, but it can also be an intracellular bacteria—proliferates within the cytosol of the host cell. The antigens of this pathogen via proteasome, endoplasmatic reticulum and Golgi apparatus are presented at the surface of this cell by the major histocompatibility complex I (MHC-I). If such a cell happens to be a dendritic cell (DC), which is an antigen-presenting cell (APC) that transports the antigen from its entrance site to the corresponding secondary lymph organ, the antigen presented can be recognized specifically by the antigen receptor (CD8) of a “matured T lymphocyte” that entered the lymphatic system. Before naive T lymphocyte have undergone maturation: first, a naïve T lymphocyte in bone marrow or thymus undergoes T-cell receptor rearrangement ( selection). T cells with high affinity to self-peptides MHC are eliminated (negative selection), whereas T cells with T-cell receptors that are able to bind self-peptides MHC molecules with at least a weak affinity survive (positive selection) and circulate in the peripheral lymphatic system. The matured T lymphocyte, recognizing the antigen by high affinity to the antigen-loaded MHC, transforms into an effector cell and proliferates. These cells are short-lived and some participate in forming memory cells. The cytotoxic T-lymphocytes (CTL) then only kill those cells, which harbor the pathogen by recognizing its antigens presented at the surface of the infected cell by MHC-I molecules. Thus, further proliferation of the pathogen is diminished. Viruses are intracellular parasites that depend on the host cell to survive and replicate. The host cell can be damaged either directly by the virus or by the immune response it provoked consisting of cytokines, macrophages, and antibodies and, most important, the CTLs. The balance of good or bad harm depends on the virus lethality, the amount of virus present (virus load), the amount of tissue infected (cyto-pathogenicity) and the affinity of CTL-response, and duration of CTL response (chronicity of the infection [10]).

One can note another path, the “exogenous-humoral” path whereby the pathogen, which is usually a bacterium, but it can be a virus or a parasite as well, proliferates in the extra-cellular space of the host. The pathogen, or fragments of it, is endocytosed into the phagolysosome of a host’s APC, which transports the antigen as a DC to the secondary lymph organ, and the antigens of the pathogen are presented at the surface of this cell by MHC-II molecules. The antigen presented can be recognized specifically by the antigen receptor of a matured helper T lymphocyte (called CD4 Th1 and CD4 Th2, resp.). A matured B lymphocyte (interacting specifically with the matured helper T lymphocyte) becomes activated (transforms into an effector cell and proliferates: these cells are short-lived, and some participate in forming memory cells), it is then called PC- (plasma-cell) producing antigen receptors (called IgG and IgE, resp.) which are soluble. These antibodies, or immune globulins, mark the pathogen, which in turn is phagozytosed and killed by macrophages.

For the function of specific adaptive/acquired immune system, the B-cell and T-cell memory is essential [11–14]. The immune memory renders the immune response at multiple encounters with the same pathogen more efficient than at the first encounter.

The mathematical model in Section 3 considers only the effector cell properties (i.e., proliferation, cell death and memory cell formation) of the immune system and the pathogen properties (i.e., proliferation, cell death and variation) thus justifying the applicability of same conceptual frame of a Darwinian entity.

##### 1.3. Previous Mathematical Models of the Immune Response

Previous approaches on the theoretical understanding of the interaction between an invading pathogen and the host’s immune system [15–23], especially on the issue of multistrain pathogens [24–27], are derived from deterministic models and are continuous in time. Continuous models provide a good representation of the dynamics when there are many participants and when fluctuations are small. These models are based on establishing a reasonable set of first-order differential equations that are assumed to be generic equations describing the properties of single cells [20]. The rates of change with respect to time of each variable describing the mean values of fractions of a total cell population are equal to a corresponding source (replication rate) and sink (death rate and rate at which new strains are generated). One studies, respectively, analytic and numerical solutions, which have mainly nonlinear properties. In the most suitable example [3, 4], these authors introduce the following differential equations with five variables as follows: where, represents the uninfected host cell, which proliferates (rate ), dies (rate ), and gets infected (rate ), represents the infected host cell, which has been infected (rate ) and dies (rates and ), represents the free virus, which proliferates within infected host cell followed by expulsion (rate ) and declines (rate ), represents the immune precursor/memory-cell, which proliferates (rate ), differentiates into immune effector cell upon antigenic challenge (rate ), and dies (rate ), represents the immune effector cell, which has differentiated from immune precursor/memory-cell (rate ) and dies (rate ).

The authors [3, 4] give parameter regions of their model, for example, the case of low virus load, where the immune system is nonresponsive, the case of high load of noncytopathic virus, where exhaustion of the immune system occurs, and the case of immune memory function where the immune response is persistent. They apply the model successfully to infections with the Lymphocyte Choriomeningitis virus (LCMV) and the HIV. Another Ansatz related to antigenic variation is given by [5] where, represents the virus variants with sequence in epitope and sequence in epitope , both coexistent, which proliferate (rate ) and being killed by CTLs (rates and ), represents the CTLs against sequence of epitope , which proliferate upon activation or being already active (rates or ) and die (rate ), represents the CTL against sequence of epitope , which proliferate upon activation or being already active (rates or ) and die (rate ).

These coupled nonlinear differential equations investigate the complex phenomena occurring in a host which is infected by a heterogeneous pathogen population, namely, inducing a fluctuating immune response against multiple epitopes with the potential of a shift of immunodominance by escape in one epitope (for a simple case the options are termed , , , and with sequences and at epitope 1 and sequences and at epitope 2, resp.).

Systems (which would die according to their differential equations approximation), when taking into account the discrete character of their microscopic components, display the emergence of macroscopic localized subpopulations with collective adaptive properties that allow their survival and development [28–30]. Simulations based on a hybrid model generate a more faithful approximation of the reality of the immune system [31].

#### 2. Developing the Methods

##### 2.1. Modeling a Darwinian Entity

Within the schema of general evolutionary biology, an entity, and thus its clonal population of individuals, undergoes a step-by-step Darwinian process from one generation to the next, that is, multiplication with random mutations and selection biased by fitness in the dependency to the actual environment (Figure 1). Each entity carries an information storage device (genotype), for example, a polymer (i.e., DNA or RNA) with a specific monomer sequence, which, in the multiplication phase, is copied with occasional mismatches (copying error probability per monomer). In the selection phase, each individual entity has a certain probability to be selected to survive according to the fitness of the phenotype (retrieved from the information storage device) in reference to its environment. Many and sustained step-by-step Darwinian processes are required from the first replicating molecule up to the emergence of mankind and many species emerged and others became extinct along the long way called Darwinian evolution.

Biological conduct is immanently stochastic, especially in the view of a cell population dynamics following a step-by-step Darwinian process. Stochastic models offer the benefit of handling the dynamics of whole population distributions (with their mean and standard deviation as deduction). These models provide a good representation of the dynamics when the numbers of participants in the process are small or when fluctuations are large. (e.g., extinction or initiation of infection). It is also worth noting that for studying extinction probabilities, it is natural to turn to stochastic models.

Some stochastic approaches deal with birth-death processes by solving “Master equations” [32], by discrete-time multitype branching processes [33, 34], and by modeling gene-amplification process with branching random walks [35, 36]. Our approach in this paper is based on the theory of branching processes, more precisely on some multitype modifications of the standard Galton-Watson processes examined in detail [37–41].

##### 2.2. Dynamical Stochastic Process of an Entity with Multiplication and Selection

This paper does not explicitly consider the information carrier (genotype) with its readouts (phenotype), nor the environment (bone marrow or thymus or secondary lymph organs in case of the immune cells, intracellular or extracellular space in case of the pathogen). The frequencies of division, the probability of forming new strains during multiplication, and the death rate, all constitute parameters in implicitly dealing with those properties. The discrete time step of the dynamics is given by the duration of each generation.

We model the process of multiplication and selection by a dynamical stochastic process with the following rules [42]: (0) Start with one individual (probability 1 of finding one individual at the end of generation ). Increase generation number from to .(i) Evaluate the probability distribution of finding individuals after multiplication phase of generation . The number of individuals reaches the cut-off value in the case of limited nutrition supply.(ii) Evaluate the probability distribution of finding individuals after selection phase of generation .(iii) Increase generation number from to and continue with (i) accordingly.

##### 2.3. Multiplication without Mutation

The probability distribution to find individuals after multiplication phase of the th generation is given by the sum over all possible paths of the conditional probabilities leading to that state ( individuals) given the state ( individuals) at the end of the selection phase of the ()th generation times the probability of that state, that is, the convolution of a binomial distribution (Figure 2(a)) where is the probability of one copy, and is the probability of no copy. The binomial coefficient counts without regard to order the number of ways of choosing copies from maximal possible copies, where is the total number of individuals after the multiplication process, and is the total number of maximal possible individuals after the multiplication process (multiplication factor ) considering the cut-off condition when the limit of nutrition supply is reached, and . is the probability of finding the population of individuals before the multiplication process (which is the same as after the selection process of ()th generation). is a probability distribution with .

The probability distribution to find individuals after the selection phase of the th generation is again given by the sum over all possible paths of the conditional probabilities leading to that state ( individuals), given the state ( individuals) at the end of the multiplication phase of the th generation times the probability of that state, that is, the convolution of a binomial distribution (Figure 2(b)) where is the probability that an individual survives, and is the probability that an individual does not survive. The binomial coefficient counts, without regard to order, the number of ways of choosing surviving individuals from a population of individuals. The probability of finding this population of individuals before the selection process is . Again, is a probability distribution with .

We assume that one initial ancestor appears in the beginning of the first generation (which is the same as after the selection process of generation ) thus the probability (root of iterative process). The numerical evaluations are given in Figure 3 (discrete limit) and Figure 4 (continuous limit). How the model responds to variations in key parameters is presented in [42]. One can insert (3) into (4) (renaming function and index, see (5)) and construct the probability generating function (PGF) with the dummy variable , see (6)

A step-by-step Darwinian process alternating between a multiplication and a selection phase per generation with a constant generation time (typically a fraction of an hour to days) and nonoverlapping generations can be contrived in the two following ways:(i)By the discrete limit with great oscillations in average numbers of individuals and (Figure 3), where each individual is copied (copy probability ), and where the selection is intermediate (selection probability to survive ). This would imply an orchestration of phases by an environmental pacemaker (day and night in case of origin of life [42–46]).(ii)By the continuous limit (Figure 4), where the probability of an individual to copy is sufficiently small (), and where the probability to be selected to survive is close to one (). This implies individual pacemakers and a smother course of the average numbers of individuals and .

Taking average values of the probability distributions, one compares the stochastic model with the deterministic model given by the differential equation [47] describing a birth-death process (number of individuals , and saturation , where ) and its analytic solution (initial value of ) Note: the probability distribution and the probability of extinction are not represented within deterministic models.

##### 2.4. Multiplication with Copying Errors Leading to One Mutant

The probability distribution to find individuals of the initial form and individuals of the mutant after the multiplication phase of the th generation is given by the convolution (Figure 5(a)) where the two remaining indices are given by Each individual (and ) replicates by a copy probability (and , resp.). (and ) are the probabilities that an error in copying the initial form occurs, and the new form emerges (and vice versa). The second binomial coefficient (and the fourth binomial coefficient) counts without regard to order the number of ways of choosing (and , resp.) error-containing copies of the initial form giving the new form (and vice versa) from a collection of total copies of the initial form (and of the new form , resp.). The first binomial coefficient (and the third binomial coefficient) counts without regard to order the number of ways of choosing total copies of the initial form (of the new form , resp.) from a collection of maximal possible copies according to the multiplication factors of the initial form (and of the new form , resp.). The cut-off conditions (where is for the total numbers of both forms and and is defined in Figure 6) then are

The probability distribution to find individuals of the initial form and individuals of the new form (mutant ) () after the selection phase of the th generation is given by the convolution (Figure 5(b)) where (and ) are the probabilities that one individual of the initial form (and mutant , resp.) survives. Figure 7 shows the dynamics resulting from evaluation by computer.

Taking average values of the probability distributions, one compares the stochastic model with the deterministic model given by the two coupled differential equations [48] being the extension of a birth-and-death process described in (7) for two entities transforming one entity into the other (number of individuals and , parameters and , and saturation ) as follows:

#### 3. Applying the Methods

##### 3.1. Modeling the Immune Response to a Pathogen by Coupling Two Darwinian Entities

Let us look first at the conceptual fundamentals of our model: to get the general idea, one applies Occam’s razor (or lex parsimoniae translating to law of succinctness) onto the complex immunological system described above and then one provides a minimal representation of the immune response. In the following we consider:(i)the host as being unstructured by not considering its multicompartmentness (i.e., not considering that the entrance site of the pathogen is spatially apart from the corresponding secondary lymph organ, where part of the immune-system response takes place);(ii)the invading pathogen () taking into account its variable antigens but not distinguishing between endogenous or exogenous paths (i.e., not considering that the pathogen thrives within a cell of the host or within the interstitial space);(iii)the immune system with the immune effector () taking into account an immune memory (), but not distinguishing between T or B lymphocytes;(iv)the step-by-step Darwinian process as fundamental to both entities (pathogen as well as to the immune effector and its memory state), which are specifically coupled;(v)the stochastic representation of a Darwinian entity as a sufficiently good starting point to solve the proposed problem.

Stochastic models offer the benefit of handling the dynamics of whole distributions with their mean and standard deviation as deduction, whereas deterministic models deal with quantities that arise as large population rescaling.

We propose a dynamical model of two interlocked Darwinian entities, the pathogen on the one hand, and the immune system on the other hand consisting of the immune effector and the immune memory (Figure 8). The coupling is such that at each time step the parameters for the pathogen system are dependent on the current state of the immune system and the parameters for the immune system are dependent on the current state of the pathogen system. As in any control system (such as body temperature of endotherms or glucose concentration in blood) there are two states: (i) the measured variable goes below a threshold or “lower set point,” then the actuator is turned on and subsequently the measured value increases, and (ii) the measured variable goes above a threshold or “upper set point,” then the actuator is turned off and subsequently the measured value decreases. Within stochastic fluctuations such systems are intrinsic periodic around a steady state.

##### 3.2. Modeling the Immune Response to One Pathogen

The multiplication phase of pathogen (with antigen ) is described by (see (3) from Section, now with index )

The selection phase of pathogen is described by (see (4) from Section, now with index )

The multiplication phase of immune effector (producing antibody ) and memory is described by (see (9) and (10) from Section 2.4, now with indices and ) where the two remaining indices are given by The cut-off conditions are

The selection phase of immune effector and memory is described by (see (12) from Section 2.4, now with indices and )

Each individual has a probability , that is, , , and , being selected to survive, it replicates by a copy-probability , multiplication-factor , that is, , , and with a mutation probability , that is, and . Thus, the multiplication factor , the error probability , and the surviving probability are parameters in function of averages and (omitting the indices and for simplicity) that form the coupling (step functions with threshold values, see captions of Figures 8, 2 and 4).

In Figure 9 the case is shown, where the pathogen is eliminated (probability of extinction after initial infection and after reinfection). In Figure 10, an oscillatory chronic case is shown, where after an apparent conquest and the subsequent relaxation of the immune reaction, the pathogen is flaring up again (probability of extinction persists below 1.0).

##### 3.3. Modeling the Immune Response to a Variable Pathogen

We consider a simple case (Figure 11) of a pathogen with two alleles at an -to- genlocus (one epitope): the pathogen (or pathogen , resp.) expressing antigen (or antigen ). The multiplication phase of pathogens and is described by ((9), (10) from Section 2.4, now with indices and ) where the two remaining indices are given by

The selection phase of pathogens and is described by ((12) from section 2.4, now with indices and )

The immune effector (or immune effector ) responding specifically, thus producing antigen-receptor a (or antigen-receptor ) which recognizes the antigen (or antigen ) and eliminates pathogen (or pathogen , resp.). While the antigen “” is the “lock” and the antigen-receptor “” is the corresponding “key” (or the antigen is the lock and the antigen-receptor is another, but corresponding key.

In Figure 12 we show a computer result of a varying pathogen ( with antigen changes into with antigen with a certain probability vice versa, see equations (20)–(22) with indices und describing pathogens and ) coupled through average values to an immune system against antigen (and against antigen ) consisting of an effector and memory (and an effector and memory ). There is no “ to ” or “ to ”—transition within the immune system. The pathogen expressing antigen is nearly eradicated, but the mutant pathogen strain-expressing antigen has escaped the immune attack (probability distribution upper left of Figure 12). As an outcome, one can see that the thriving of the pathogen within the host is prolonged.

#### 4. Maturation of T-Lymphocytes

As mentioned in Section 1.2, the T-lymphocyte comes in four different forms: a naïve T-lymphocyte in bone marrow or thymus undergoes T-cell receptor rearrangement (-selection), where T-cells with high affinity to self-peptides MHC are eliminated (negative selection), and T-cells with T cell receptors that are able to bind MHC molecules with at least a weak affinity survive in the peripheral lymphatic system (positive selection). The matured T-lymphocyte recognizing the antigen by high affinity to the antigen loaded MHC transforms into an effector cell and proliferates. We consider in Figure 13 the dynamics of a probability function with four variables describing naïve lymphocyte , lymphocyte with strong affinity to self-peptides, matured lymphocyte with weak affinity to foreign-peptides, this lymphocyte gets inactivated, matured lymphocyte with strong affinity to foreign-peptides, this lymphocyte gets activated (Mat = matured, = inactivated, = activated). The multiplication phase is described by where the remaining indices are given by For the phase () in Figure 13, the multiplication phase is described by inserting formula (25) into formula (23) where the remaining indices instead of (24) are given by