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

Authors: Christoph Kuhn
Journal: Computational and Mathematical Methods in Medicine
 (2012)

---

## 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.


---

## Body

## 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 (Section1.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 Section3 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:
(1)x˙=λ-(d-βv)x,y˙=βxv-(a+pz)y,v˙=ky-uv,w˙=(cy-fy-r)w,z˙=fyw-bz,
where, x represents the uninfected host cell, which proliferates (rate λ), dies (rate dx),  and gets infected (rate βvx), y represents the infected host cell, which has been infected (rate βxv) and dies (rates ay and pzy), v represents the free virus, which proliferates within infected host cell followed by expulsion (rate ky) and declines (rate uv), w represents the immune precursor/memory-cell, which proliferates (rate cyw), differentiates into immune effector cell upon antigenic challenge (rate fyw), and dies (rate rw), z represents the immune effector cell, which has differentiated from immune precursor/memory-cell (rate fyw) and dies (rate bz).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]
(2)v˙ij=(rij-pixi-qjyj)vij,x˙i=ηci∑jvij+(ci∑jvij-b)xi,y˙j=ηkj∑jvij+(kj∑jvij-b)yi,
where, vij represents the virus variants with sequence i in epitope A and sequence j in epitope B, both coexistent, which proliferate (rate rij) and being killed by CTLs (rates pixivij and qjyjvij), xi represents the CTLs against sequence i of epitope A, which proliferate  upon activation or being already active (rates ηci∑jvij or ci(∑jvij)xi) and die (rate bxi), yj represents the CTL against sequence j of epitope B, which proliferate upon activation or being already active (rates ηkj∑jvij or kj(∑jvij)yj) and die (rate byj).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 termedA1, B1, C2, and D2 with sequences A and B at epitope 1 and sequences C and D 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].


## 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 (Section1.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 Section3 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:
(1)x˙=λ-(d-βv)x,y˙=βxv-(a+pz)y,v˙=ky-uv,w˙=(cy-fy-r)w,z˙=fyw-bz,
where, x represents the uninfected host cell, which proliferates (rate λ), dies (rate dx),  and gets infected (rate βvx), y represents the infected host cell, which has been infected (rate βxv) and dies (rates ay and pzy), v represents the free virus, which proliferates within infected host cell followed by expulsion (rate ky) and declines (rate uv), w represents the immune precursor/memory-cell, which proliferates (rate cyw), differentiates into immune effector cell upon antigenic challenge (rate fyw), and dies (rate rw), z represents the immune effector cell, which has differentiated from immune precursor/memory-cell (rate fyw) and dies (rate bz).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]
(2)v˙ij=(rij-pixi-qjyj)vij,x˙i=ηci∑jvij+(ci∑jvij-b)xi,y˙j=ηkj∑jvij+(kj∑jvij-b)yi,
where, vij represents the virus variants with sequence i in epitope A and sequence j in epitope B, both coexistent, which proliferate (rate rij) and being killed by CTLs (rates pixivij and qjyjvij), xi represents the CTLs against sequence i of epitope A, which proliferate  upon activation or being already active (rates ηci∑jvij or ci(∑jvij)xi) and die (rate bxi), yj represents the CTL against sequence j of epitope B, which proliferate upon activation or being already active (rates ηkj∑jvij or kj(∑jvij)yj) and die (rate byj).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 termedA1, B1, C2, and D2 with sequences A and B at epitope 1 and sequences C and D 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 (Figure1).  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.Figure 1
Schema of a Darwinian entity. An individual is singled out from the population. A period of one generation is shown. Incidental copying error occurs during multiplication (changed genotype resulting in changed phenotype) with new fitness in reference to its structured environment. The probability of being selected to survive is given according to new fitness.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 individualW0S(1)=1 (probability 1 of finding one individual at the end of generation n=0). Increase generation number from n=0 to n=1.(i)
Evaluate the probability distributionWnM(ν) of finding 0≤ν≤Nmax individuals after multiplication phase M of generation n. The number of individuals reaches the cut-off value Nmax in the case of limited nutrition supply.(ii)
Evaluate the probability distributionWnS(N) of finding 0≤N≤Nmax individuals after selection phase S of generation n.(iii)
Increase generation number fromn to n+1 and continue with (i) accordingly.


### 2.3. Multiplication without Mutation

The probability distribution to find0≤ν≤Nmax individuals after multiplication phase M of the nth generation is given by the sum over all possible paths of the conditional probabilities (μ-Nν-N)α(ν-N)(1-α)(μ-ν) leading to that state (ν individuals) given the state (N individuals) at the end of the selection phase of the (n-1)th generation times the probability Wn-1S(N) of that state, that is, the convolution of a binomial distribution (Figure 2(a))
(3)WnM(ν)=∑N=ην(μ-Nν-N)α(ν-N)(1-α)(μ-ν)·Wn-1S(N),
where α is the probability of one copy, and 1-α is the probability of no copy. The binomial coefficient counts without regard to order the number of ways of choosing ν-N copies from μ-N maximal possible copies, where ν is the total number of individuals after the multiplication process, and μ=Min(ρ·N,Nmax) 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 η=Ceiling(ν/ρ). Wn-1S(N) is the probability of finding the population of N individuals before the multiplication process (which is the same as after the selection process of (n-1)th generation). WnM(ν) is a probability distribution with ∑ν=0NmaxWnM(ν)=1.Figure 2
Sketch of how convolution of a binomial distribution is applied to probability distribution. (a) multiplication phaseM, (3). (b) Selection phase S, (4). (μ-N) Maximal possible number of copies. (ν-N) Number of copies.The probability distribution to find0≤N≤Nmax individuals after the selection phase S of the nth generation is again given by the sum over all possible paths of the conditional probabilities (νN)βN(1-β)(ν-N) leading to that state (N individuals), given the state (ν individuals) at the end of the multiplication phase of the nth generation times the probability WnM(ν) of that state, that is, the convolution of a binomial distribution (Figure 2(b))
(4)WnS(N)=∑ν=0Nmax(νN)βN(1-β)(ν-N)·WnM(ν),
where β is the probability that an individual survives, and 1-β is the probability that an individual does not survive. The binomial coefficient counts, without regard to order, the number of ways of choosing N surviving individuals from a population of ν individuals. The probability of finding this population of n individuals before the selection process is WnM(ν). Again, WnS(N) is a probability distribution with ∑N=0NmaxWnS(N)=1.We assume that one initial ancestor appears in the beginning of the first generation (which is the same as after the selection process of generationn=0) thus the probability W0S(1)=1 (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 s, see (6)
(5)PN(j)=∑ν=jμ(μ-Nν-N)α(ν-N)(1-α)(μ-ν)·(νj)βj(1-β)(ν-j),(6)GN(s)=∑j=0NmaxPN(j)·sj=(1-β+βs)N(1-αβ+αβs)N.Figure 3
Dynamical stochastic process of multiplication and selection. Discrete limit, one initial ancestor. The two cardinal examples                       (a)n=7 (rn-regime, r=ρβ=4/3, ln(r)=0.288) and (b) n=50 (Z∞-regime). Upper left: probability distribution WnM(ν) of finding 0≤ν≤Nmax individuals after the multiplication phase of the nth generation (3). Lower left probability distribution WnS(N) of finding 0≤N≤Nmax individuals after the selection phase of the nth generation (4). Maximal total number Nmax (cut-off value due to limited nutrition supply). Upper right: probability of extinction W(0) along generation n. Lower right: average number of individuals ν-=∑ν=0NmaxWnM(ν)·ν and N-=∑N=0NmaxWnS(N)·N along generation n. Parameters: Nmax=128, multiplication factor ρ=2, and copy probability α=1 (all individuals that are present before multiplication copy once), surviving probability β=2/3. Deterministic model (green in lower right, equations (7) and (8)): K=64, N0=1, and R=0.344.Figure 4
Dynamical stochastic process of multiplication and selection. Continuous limit, one initial ancestor. Two cardinal examples (a)n=33 and (b) n=150. Upper left: probability distribution WnM(ν) of finding ν individuals after multiplication phase M of generation n (3). Lower left: probability distribution WnS(N) of finding N individuals after selection phase S of generation n (4). Average and standard deviation indicated (correct values only by taking the probability of extinction W(0) into account). Upper right: probability of extinction W(0) along generation n. Lower right: number of individuals ν- and N- with its standard deviation along generation n (least square fit, R=0.0694, see (7) and (8)). Parameters: Maximal total number  Nmax=128, multiplication factor ρ=2, and copy probability α=0.1 (only 10% of individuals that are present before multiplication copy once), surviving probability β=0.9664.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 N- (Figure 3), where each individual is copied (copy probability α=1), and where the selection is intermediate (selection probability to survive β=2/3). 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 (Figure4), where the probability of an individual to copy is sufficiently small (α=0.1), and where the probability to be selected to survive is close to one (β=0.9664). This implies individual pacemakers and a smother course of the average numbers of individuals ν- and N-.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 N, R=birthrate-deathrate and saturation K, where dN/dt=0) and its analytic solution (initial value N0 of N)
(7)dNdt=RN-RKN2,(8)N(t)=eRtN0K(K-N0+N0eRt).
Note: the probability distribution W(N) and the probability of extinction W(0) are not represented within deterministic models.


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

The probability distribution to findνA individuals of the initial form A and νB individuals of the mutant B(0≤νA+νB=ν≤Nmax) after the multiplication phase M of the nth generation is given by the convolution (Figure 5(a))
(9)WnM(νA,νB)=∑NA=0Nmax∑NB=0Nmax-NA∑kA=0μA-NA∑kAB=0μA-kA-NAWn-1S(NA,NB)·(μA-NAkA)αA(μA-kA-NA)(1-αA)kA·(μA-kA-NAkAB)εABkAB(1-εAB)(μA-kA-kAB-NA)·(μB-NBkB)αB(μB-kB-NB)(1-αB)kB·(μB-kB-NBkBA)εBAkBA(1-εBA)(μB-kB-kBA-NB),
where the two remaining indices are given by
(10)kB=μA+μB-νA-νB-kA,kBA=νA+kAB+kA-μA.
Each individual A (and B) replicates by a copy probability αA (and αB, resp.). εAB (and εBA) are the probabilities that an error in copying the initial form A occurs, and the new form B emerges (and vice versa). The second binomial coefficient (and the fourth binomial coefficient) counts without regard to order the number of ways of choosing kAB (and kBA, resp.) error-containing copies of the initial form A giving the new form B (and vice versa) from a collection of  μA-kA-NA(andμB-kB-NB) total copies of the initial form A (and of the new form B, resp.). The first binomial coefficient (and the third binomial coefficient) counts without regard to order the number of ways of choosing μA-kA-NA(andμB-kB-NB) total copies of the initial form A (of the new form B, resp.) from a collection of μA-NA(andμB-NB) maximal possible copies according to the multiplication factors ρA(andρB) of the initial form A (and of the new form B, resp.). The cut-off conditions (where Nmax is for the total numbers of both forms A and B and νadj is defined in Figure 6) then are
(11)μA=ρA·NAμB=ρB·NBifρA·NA+ρB·NB≤Nmax,μA=νadjμB=Nmax-νadjifρA·NA+ρB·NB>Nmax.Figure 5
Sketch for (9), (10), and (12) of how convolution of a binomial distribution is applied to probability distribution.                      (a) Multiplication phase M, (9), (10): μA-NA-kA total copies of initial form A (and μB-NB-kB total copies of mutant B); μA-NA maximal possible copies according to the multiplication factors ρA of initial form A (and μB-NB maximal possible copies according to the multiplication factors ρB of mutant B); kAB (and kBA) number of copies transforming from initial form A to mutant B (and vice versa resp.).  (b) Selection phase S, (12).Figure 6
Sketch to construct the cut-off condition in the case of initial formA and one kind of mutant B. If, by multiplication of {NA,NB} individuals (black disk •) with factors ρA and ρB respectively, the total number ρA·NA+ρB·NB>Nmax would be beyond the limit of supply Nmax, (circle ∘), the correct cut-off point {νadj,Nmax-νadj} (black square ■) is the most adjacent to the intersection. Note: partition and take average for more than one such point. The probability distributions and cut-off conditions in the cases of more than one kind of mutant are accordingly.The probability distribution to findNA individuals of the initial form A and NB individuals of the new form (mutant B) (0≤NA+NB=N≤Nmax) after the selection phase S of the nth generation is given by the convolution (Figure 5(b))
(12)WnS(NA,NB)=∑νA=NANmax∑νB=NBNmax-νAWnM(νA,νB)·(νANA)βANA(1-βA)(νA-NA)·(νBNB)βBNB(1-βB)(νB-NB),
where βA (and βB) are the probabilities that one individual of the initial form A (and mutant B, resp.) survives. Figure 7 shows the dynamics resulting from evaluation by computer.Figure 7
Dynamical stochastic process for initial formA and mutant B. One initial ancestor of form A. Three cardinal examples (a) n=48, (b) n=164, and (c) n=600. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for n): probability distribution WnM(νA,νB) of finding νA individuals of initial form A (red) and νB individuals of mutant B (blue) after multiplication phase M of generation n (9), (10). Lower left (triangle graph): probability distribution WnS(NA,NB) of finding NA individuals of initial form A (red) and NB individuals of mutant B (blue) after selection phase S of generation n (12). Upper right: total extinction probability W(0) of both initial form (A) and mutant (B) together (violet) along generation n. Lower right: average number of individuals ν- and N- of initial form A (red) and mutant B (blue) with their standard deviations (black) along generation n. Parameters: maximal total number Nmax=32, copy probability αA=0.1 and αB=0.1, multiplication factor ρA=3 of initial form A, multiplication factor ρB=6 of mutant B, mutation probability εAB=0.01 from initial form A to mutant B, mutation probability εBA=0.001 from mutant B to initial form A, surviving probability βA=0.95 of initial form A, surviving probability βB=0.95 of mutant B. Deterministic model (green in lower right, (13)): Ra=0.128, Rb=0.14, and K=16.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 NA and NB, R=birthrate-deathrate parameters Ra and Rb, and saturation K) as follows:
(13)dNAdt=Ra·NA-RaKNA2-RbKNANB,dNBdt=Rb·NB-RbKNB2-RaKNBNA.


## 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 (Figure1).  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.Figure 1
Schema of a Darwinian entity. An individual is singled out from the population. A period of one generation is shown. Incidental copying error occurs during multiplication (changed genotype resulting in changed phenotype) with new fitness in reference to its structured environment. The probability of being selected to survive is given according to new fitness.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 individualW0S(1)=1 (probability 1 of finding one individual at the end of generation n=0). Increase generation number from n=0 to n=1.(i)
Evaluate the probability distributionWnM(ν) of finding 0≤ν≤Nmax individuals after multiplication phase M of generation n. The number of individuals reaches the cut-off value Nmax in the case of limited nutrition supply.(ii)
Evaluate the probability distributionWnS(N) of finding 0≤N≤Nmax individuals after selection phase S of generation n.(iii)
Increase generation number fromn to n+1 and continue with (i) accordingly.


## 2.3. Multiplication without Mutation

The probability distribution to find0≤ν≤Nmax individuals after multiplication phase M of the nth generation is given by the sum over all possible paths of the conditional probabilities (μ-Nν-N)α(ν-N)(1-α)(μ-ν) leading to that state (ν individuals) given the state (N individuals) at the end of the selection phase of the (n-1)th generation times the probability Wn-1S(N) of that state, that is, the convolution of a binomial distribution (Figure 2(a))
(3)WnM(ν)=∑N=ην(μ-Nν-N)α(ν-N)(1-α)(μ-ν)·Wn-1S(N),
where α is the probability of one copy, and 1-α is the probability of no copy. The binomial coefficient counts without regard to order the number of ways of choosing ν-N copies from μ-N maximal possible copies, where ν is the total number of individuals after the multiplication process, and μ=Min(ρ·N,Nmax) 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 η=Ceiling(ν/ρ). Wn-1S(N) is the probability of finding the population of N individuals before the multiplication process (which is the same as after the selection process of (n-1)th generation). WnM(ν) is a probability distribution with ∑ν=0NmaxWnM(ν)=1.Figure 2
Sketch of how convolution of a binomial distribution is applied to probability distribution. (a) multiplication phaseM, (3). (b) Selection phase S, (4). (μ-N) Maximal possible number of copies. (ν-N) Number of copies.The probability distribution to find0≤N≤Nmax individuals after the selection phase S of the nth generation is again given by the sum over all possible paths of the conditional probabilities (νN)βN(1-β)(ν-N) leading to that state (N individuals), given the state (ν individuals) at the end of the multiplication phase of the nth generation times the probability WnM(ν) of that state, that is, the convolution of a binomial distribution (Figure 2(b))
(4)WnS(N)=∑ν=0Nmax(νN)βN(1-β)(ν-N)·WnM(ν),
where β is the probability that an individual survives, and 1-β is the probability that an individual does not survive. The binomial coefficient counts, without regard to order, the number of ways of choosing N surviving individuals from a population of ν individuals. The probability of finding this population of n individuals before the selection process is WnM(ν). Again, WnS(N) is a probability distribution with ∑N=0NmaxWnS(N)=1.We assume that one initial ancestor appears in the beginning of the first generation (which is the same as after the selection process of generationn=0) thus the probability W0S(1)=1 (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 s, see (6)
(5)PN(j)=∑ν=jμ(μ-Nν-N)α(ν-N)(1-α)(μ-ν)·(νj)βj(1-β)(ν-j),(6)GN(s)=∑j=0NmaxPN(j)·sj=(1-β+βs)N(1-αβ+αβs)N.Figure 3
Dynamical stochastic process of multiplication and selection. Discrete limit, one initial ancestor. The two cardinal examples                       (a)n=7 (rn-regime, r=ρβ=4/3, ln(r)=0.288) and (b) n=50 (Z∞-regime). Upper left: probability distribution WnM(ν) of finding 0≤ν≤Nmax individuals after the multiplication phase of the nth generation (3). Lower left probability distribution WnS(N) of finding 0≤N≤Nmax individuals after the selection phase of the nth generation (4). Maximal total number Nmax (cut-off value due to limited nutrition supply). Upper right: probability of extinction W(0) along generation n. Lower right: average number of individuals ν-=∑ν=0NmaxWnM(ν)·ν and N-=∑N=0NmaxWnS(N)·N along generation n. Parameters: Nmax=128, multiplication factor ρ=2, and copy probability α=1 (all individuals that are present before multiplication copy once), surviving probability β=2/3. Deterministic model (green in lower right, equations (7) and (8)): K=64, N0=1, and R=0.344.Figure 4
Dynamical stochastic process of multiplication and selection. Continuous limit, one initial ancestor. Two cardinal examples (a)n=33 and (b) n=150. Upper left: probability distribution WnM(ν) of finding ν individuals after multiplication phase M of generation n (3). Lower left: probability distribution WnS(N) of finding N individuals after selection phase S of generation n (4). Average and standard deviation indicated (correct values only by taking the probability of extinction W(0) into account). Upper right: probability of extinction W(0) along generation n. Lower right: number of individuals ν- and N- with its standard deviation along generation n (least square fit, R=0.0694, see (7) and (8)). Parameters: Maximal total number  Nmax=128, multiplication factor ρ=2, and copy probability α=0.1 (only 10% of individuals that are present before multiplication copy once), surviving probability β=0.9664.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 N- (Figure 3), where each individual is copied (copy probability α=1), and where the selection is intermediate (selection probability to survive β=2/3). 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 (Figure4), where the probability of an individual to copy is sufficiently small (α=0.1), and where the probability to be selected to survive is close to one (β=0.9664). This implies individual pacemakers and a smother course of the average numbers of individuals ν- and N-.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 N, R=birthrate-deathrate and saturation K, where dN/dt=0) and its analytic solution (initial value N0 of N)
(7)dNdt=RN-RKN2,(8)N(t)=eRtN0K(K-N0+N0eRt).
Note: the probability distribution W(N) and the probability of extinction W(0) are not represented within deterministic models.


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

The probability distribution to findνA individuals of the initial form A and νB individuals of the mutant B(0≤νA+νB=ν≤Nmax) after the multiplication phase M of the nth generation is given by the convolution (Figure 5(a))
(9)WnM(νA,νB)=∑NA=0Nmax∑NB=0Nmax-NA∑kA=0μA-NA∑kAB=0μA-kA-NAWn-1S(NA,NB)·(μA-NAkA)αA(μA-kA-NA)(1-αA)kA·(μA-kA-NAkAB)εABkAB(1-εAB)(μA-kA-kAB-NA)·(μB-NBkB)αB(μB-kB-NB)(1-αB)kB·(μB-kB-NBkBA)εBAkBA(1-εBA)(μB-kB-kBA-NB),
where the two remaining indices are given by
(10)kB=μA+μB-νA-νB-kA,kBA=νA+kAB+kA-μA.
Each individual A (and B) replicates by a copy probability αA (and αB, resp.). εAB (and εBA) are the probabilities that an error in copying the initial form A occurs, and the new form B emerges (and vice versa). The second binomial coefficient (and the fourth binomial coefficient) counts without regard to order the number of ways of choosing kAB (and kBA, resp.) error-containing copies of the initial form A giving the new form B (and vice versa) from a collection of  μA-kA-NA(andμB-kB-NB) total copies of the initial form A (and of the new form B, resp.). The first binomial coefficient (and the third binomial coefficient) counts without regard to order the number of ways of choosing μA-kA-NA(andμB-kB-NB) total copies of the initial form A (of the new form B, resp.) from a collection of μA-NA(andμB-NB) maximal possible copies according to the multiplication factors ρA(andρB) of the initial form A (and of the new form B, resp.). The cut-off conditions (where Nmax is for the total numbers of both forms A and B and νadj is defined in Figure 6) then are
(11)μA=ρA·NAμB=ρB·NBifρA·NA+ρB·NB≤Nmax,μA=νadjμB=Nmax-νadjifρA·NA+ρB·NB>Nmax.Figure 5
Sketch for (9), (10), and (12) of how convolution of a binomial distribution is applied to probability distribution.                      (a) Multiplication phase M, (9), (10): μA-NA-kA total copies of initial form A (and μB-NB-kB total copies of mutant B); μA-NA maximal possible copies according to the multiplication factors ρA of initial form A (and μB-NB maximal possible copies according to the multiplication factors ρB of mutant B); kAB (and kBA) number of copies transforming from initial form A to mutant B (and vice versa resp.).  (b) Selection phase S, (12).Figure 6
Sketch to construct the cut-off condition in the case of initial formA and one kind of mutant B. If, by multiplication of {NA,NB} individuals (black disk •) with factors ρA and ρB respectively, the total number ρA·NA+ρB·NB>Nmax would be beyond the limit of supply Nmax, (circle ∘), the correct cut-off point {νadj,Nmax-νadj} (black square ■) is the most adjacent to the intersection. Note: partition and take average for more than one such point. The probability distributions and cut-off conditions in the cases of more than one kind of mutant are accordingly.The probability distribution to findNA individuals of the initial form A and NB individuals of the new form (mutant B) (0≤NA+NB=N≤Nmax) after the selection phase S of the nth generation is given by the convolution (Figure 5(b))
(12)WnS(NA,NB)=∑νA=NANmax∑νB=NBNmax-νAWnM(νA,νB)·(νANA)βANA(1-βA)(νA-NA)·(νBNB)βBNB(1-βB)(νB-NB),
where βA (and βB) are the probabilities that one individual of the initial form A (and mutant B, resp.) survives. Figure 7 shows the dynamics resulting from evaluation by computer.Figure 7
Dynamical stochastic process for initial formA and mutant B. One initial ancestor of form A. Three cardinal examples (a) n=48, (b) n=164, and (c) n=600. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for n): probability distribution WnM(νA,νB) of finding νA individuals of initial form A (red) and νB individuals of mutant B (blue) after multiplication phase M of generation n (9), (10). Lower left (triangle graph): probability distribution WnS(NA,NB) of finding NA individuals of initial form A (red) and NB individuals of mutant B (blue) after selection phase S of generation n (12). Upper right: total extinction probability W(0) of both initial form (A) and mutant (B) together (violet) along generation n. Lower right: average number of individuals ν- and N- of initial form A (red) and mutant B (blue) with their standard deviations (black) along generation n. Parameters: maximal total number Nmax=32, copy probability αA=0.1 and αB=0.1, multiplication factor ρA=3 of initial form A, multiplication factor ρB=6 of mutant B, mutation probability εAB=0.01 from initial form A to mutant B, mutation probability εBA=0.001 from mutant B to initial form A, surviving probability βA=0.95 of initial form A, surviving probability βB=0.95 of mutant B. Deterministic model (green in lower right, (13)): Ra=0.128, Rb=0.14, and K=16.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 NA and NB, R=birthrate-deathrate parameters Ra and Rb, and saturation K) as follows:
(13)dNAdt=Ra·NA-RaKNA2-RbKNANB,dNBdt=Rb·NB-RbKNB2-RaKNBNA.


## 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 (P) 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 (E) taking into account an immune memory (M), 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 pathogenP on the one hand, and the immune system on the other hand consisting of the immune effector E and the immune memory M (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.Figure 8
Scheme of the immune system response (immune effectorsEa producing antigen-receptor a and immune-memory cells Ma) to an invading pathogen (pathogen PA carrying antigen A). Two-state model: immune effector is turned on (from neutralized state to effector state) when average number of pathogens go above threshold TP; pathogen is turned off (from proliferating to neutralized state) when average immune response goes above threshold TE; immune effector is turned off (from effector state to neutralized state) when average number of pathogens go below threshold. First infection by the proliferating pathogen PA (high multiplication factor, high selection probability to survive) initiates the immune response of the host: an immune cell with the recipe for antigen-receptor “a” is singled out from the reservoir of naive immune cells to form an immune effector Ea which proliferates (high-multiplication factor, high-selection probability). Above a threshold titer of antigen-receptor a, the pathogen is neutralized PA (low multiplication factor, low selection probability). The immune effectors Ea are transformed into immune-memory cells Ma which do not produce antigen-receptor a but carry its recipe (low-multiplication factor, high-selection probability). During any further infections by the pathogen PA, the immune-memory cells with the recipe for antigen-receptor a are formed back into immune effectors producing antigen-receptor a. Pathogen PA carrying antigen A can transform to pathogen PB carrying antigen B. A new immune respond has to be launched with immune effectors Eb producing antigen-receptor b and immune-memory cells Mb.


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

The multiplication phase of pathogenP (with antigen A) is described by (see (3) from Section, now with index P)
(14)WnM(νP)=∑NP=ην(μP-NPνP-NP)αP(νP-NP)(1-αP)(μP-νP)·Wn-1S(NP).The selection phase of pathogenP is described by (see (4) from Section, now with index P)
(15)WnS(NP)=∑ν=NPNmax(νPNP)βPNP(1-βP)(νP-NP)·WnM(νP).The multiplication phase of immune effectorE (producing antibody a) and memory M is described by (see (9) and (10) from Section 2.4, now with indices E and M)
(16)WnM(νE,νM)=∑NE=0Nmax∑NM=0Nmax-NE∑kE=0μE-NE∑kEM=0μE-kE-NEWn-1S(NE,NM)·(μE-NEkE)αE(μE-kE-NE)(1-αE)kE·(μE-kE-NEkEM)εEMkEM(1-εEM)(μE-kE-kEM-NE)·(μM-NMkM)αM(μM-kM-NM)(1-αM)kM·(μM-kM-NMkME)εMEkME(1-εME)(μM-kM-kME-NM),
where the two remaining indices are given by
(17)kM=μE+μM-νE-νM-kE,kME=νE+kEM+kE-μE.
The cut-off conditions are
(18)μE=ρE·NEμM=ρM·NMifρE·NE+ρM·NM≤Nmax,μE=νadjμM=Nmax-νadjifρE·NE+ρM·NM>Nmax,μP=ρP·NPifρP·NP≤Nmax,μP=NmaxifρP·NP>Nmax.The selection phase of immune effectorE and memory M is described by (see (12) from Section 2.4, now with indices E and M)
(19)WnS(NE,NM)=∑νE=NENmax∑νM=NMNmax-νEWnM(νE,νM)·(νENE)βENE(1-βE)(νE-NE)·(νMNM)βMNM(1-βM)(νM-NM).Each individual has a probabilityβ, that is, βP(E-), βE(P-), and βM(P-), being selected to survive, it replicates by a copy-probability α=0.1, multiplication-factor ρ, that is, ρP(E-), ρE(P-), and ρM(P-) with a mutation probability ε, that is, εEM(P-) and εME(P-). Thus, the multiplication factor ρ, the error probability ε, and the surviving probability β are parameters in function of averages P-=∑NP=0NP,maxW(NP)·NP and E-=∑NE,NM=0NmaxW(NE,NM)·NE (omitting the indices S/M and n for simplicity) that form the coupling (step functions with threshold values, see captions of Figures 8, 2 and 4).In Figure9 the case is shown, where the pathogen is eliminated (probability of extinction W(0)=1.0 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 W(0) persists below 1.0).Figure 9
Computer result of a pathogenP (with antigen A) coupled through averages to an immune system consisting of an effector E (producing antigen receptor a) and memory M. Case where pathogen is eliminated. Glimpses at generation (a) n=13, (b) n=84 and (c) n=200. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for N): probability WnS(NP) of finding 0≤NP≤Nmax individuals of the pathogen P (green). Lower left (triangle graph): probability WnS(NE,NM) of finding NE individuals of the immune effector E (red) and NM individuals of the immune memory M (blue), (0≤NE+NM≤Nmax). Upper right: extinction probability W(0) as a function of generations n of pathogen (green), immune effector, and immune memory (violet). Lower right: average of pathogen P- (green), of immune effector E- (red), and of immune memory M- (blue) as a function of generations n. Parameter values: maximal total number Nmax=32; αP=0.1, αE=0.1, and αM=0.1. For average E-<TE=1.5 (below threshold value, immune effector E inactive): ρP=6, βP=0.96. For average E->TE=1.5 (above threshold value, immune effector E active): ρP=2, βP=0.65. For average P->TP=0.5 (below threshold value, pathogen P not seen by immune system): ρE=2, ρM=2, εEM=0.50, εME=0.01, βE=0.90, and βM=0.93. For average P->TP=0.5 (above threshold value, pathogen P seen by immune system): ρE=6, ρM=2, εEM=0.50, εME=0.50, βE=0.95, and βM=0.93.Figure 10
Computer result of a pathogenP (with antigen A) coupled through averages to an immune system consisting of an effector E (producing antigen receptor a) and memory M. Case where the pathogen is reappearing while the immune system response is low. Glimpse at generation n=100. Parameter values other than Figure 2: For average E->TE=1.5 (above threshold value, immune effector E active): βP=0.75. For average P-<TP=1.5 (below threshold value, pathogen P not seen by immune system): βE=0.85. For average  P->TP=1.5 (above threshold value, pathogen P seen by immune system): βE=0.90.


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

We consider a simple case (Figure11) of a pathogen with two alleles at an A-to-B genlocus (one epitope): the pathogen PA (or pathogen PB, resp.) expressing antigen A (or antigen B).        The multiplication phase of pathogens PA and PB is described by ((9), (10) from Section 2.4, now with indices PA and PB)
(20)WnM(νPA,νPB)=∑NPA=0Nmax∑NPB=0Nmax-NPA∑kPA=0μPA-NPA∑kPAPB=0μPA-kPA-NPAWn-1S(NPA,NPB).(μPA-NPAkPA)αPA(μPA-kPA-NPA)(1-αPA)kPA.(μPA-kPA-NPAkPAPB)εPAPBkPAPB(1-εPAPB)(μPA-kPA-kPAPB-NPA).(μPB-NPBkPB)αPB(μPB-kPB-NPB)(1-αPB)kPB.(μPB-kPB-NPBkPBPA)εPBPAkPBPA(1-εPBPA)(μPB-kPB-kPBPA-NPB),
where the two remaining indices are given by
(21)kPB=μPA+μPB-νPA-νPB-kPA,kPBPA=νPA+kPAPB+kPA-μPA.Figure 11
Immune response to a variable pathogen (pathogen strains each with different antigens which are presented to the immune system). Two alleles at aAB-genlocus of the pathogen expressing antigen A or B: the immune effectors respond specifically (lock and key principle) by proliferation and producing antigen-receptor (antibody) a or b.The selection phase of pathogensPA and PB is described by ((12) from section 2.4, now with indices PA and PB)
(22)WnS(NPA,NPB)=∑νPA=NPANmax∑νPB=NPBNmax-νPAWnM(νPA,νPB).(νPANPA)βPANPA(1-βPA)(νPA-NPA).(νPBNPB)βPBNPB(1-βPB)(νPB-NPB).The immune effectorEa (or immune effector Eb) responding specifically, thus producing antigen-receptor a (or antigen-receptor b) which recognizes the antigen A (or antigen B) and eliminates pathogen PA (or pathogen PB, resp.). While the antigen “A” is the “lock” and the antigen-receptor “a” is the corresponding “key” (or the antigen B is the lock and the antigen-receptor b is another, but corresponding key.In Figure12 we show a computer result of a varying pathogen (PA with antigen A changes into PB with antigen B with a certain probability εPAPB,εPBPA vice versa, see equations (20)–(22) with indices PA und PB describing pathogens A and B) coupled through average values to an immune system against antigen A (and against antigen B) consisting of an effector Ea and memory Ma (and an effector Eb and memory Mb). There is no “a to b” or “b to a”—transition within the immune system. The pathogen expressing antigen A is nearly eradicated, but the mutant pathogen strain-expressing antigen B 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.Figure 12
Computer result of a varying pathogen (PA with antigen A changing into PB with antigen B) coupled through averages to an immune system against antigen A consisting of an effector Ea and memory Ma and against antigen B consisting of an effector Eb and memory Mb. Partial change of pathogen PA to pathogen PB escaping immune effector Ea and with delayed immune response of effector Eb. Glimpses at generation (a) n=22, (b) n=100. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for N): probability WnS(NPA,NPB) of finding NPA individuals of the pathogen PA with antigen A and of finding NPB individuals of the pathogen PB with antigen B (green) (0≤NPA+NPB≤Nmax). Middle left (triangle graph): probability WnS(NEa,NMa) of finding NEa individuals of the immune effector (red) and NMa individuals of immune memory (blue) against antigen A (0≤NEa+NMa≤Nmax). Lower left (triangle graph): probability WnS(NEb,NMb) of finding NEb individuals of the immune effector (red) and NMb individuals of immune memory (blue) against antigen B (0≤NEb+NMb≤Nmax). Upper right: extinction probability W(0) as a function of generations n of pathogen (green), immune effector & memory against antigen A and immune effector & memory against antigen B (both violet). Middle and lower right: average of pathogen P- (green), of immune effector E- (red) and of immune memory M- (blue) against antigen A (middle right) and against antigen B (lower right) as a function of generations n. Parameter values others than Figure 2: average TPA=1.5  and average TPB=1.5, respectively (threshold value to switch immune system); εPAPB=0.01, εPBPA=0.001 (pathogen variability).


## 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 (P) 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 (E) taking into account an immune memory (M), 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 pathogenP on the one hand, and the immune system on the other hand consisting of the immune effector E and the immune memory M (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.Figure 8
Scheme of the immune system response (immune effectorsEa producing antigen-receptor a and immune-memory cells Ma) to an invading pathogen (pathogen PA carrying antigen A). Two-state model: immune effector is turned on (from neutralized state to effector state) when average number of pathogens go above threshold TP; pathogen is turned off (from proliferating to neutralized state) when average immune response goes above threshold TE; immune effector is turned off (from effector state to neutralized state) when average number of pathogens go below threshold. First infection by the proliferating pathogen PA (high multiplication factor, high selection probability to survive) initiates the immune response of the host: an immune cell with the recipe for antigen-receptor “a” is singled out from the reservoir of naive immune cells to form an immune effector Ea which proliferates (high-multiplication factor, high-selection probability). Above a threshold titer of antigen-receptor a, the pathogen is neutralized PA (low multiplication factor, low selection probability). The immune effectors Ea are transformed into immune-memory cells Ma which do not produce antigen-receptor a but carry its recipe (low-multiplication factor, high-selection probability). During any further infections by the pathogen PA, the immune-memory cells with the recipe for antigen-receptor a are formed back into immune effectors producing antigen-receptor a. Pathogen PA carrying antigen A can transform to pathogen PB carrying antigen B. A new immune respond has to be launched with immune effectors Eb producing antigen-receptor b and immune-memory cells Mb.


## 3.2. Modeling the Immune Response to One Pathogen

The multiplication phase of pathogenP (with antigen A) is described by (see (3) from Section, now with index P)
(14)WnM(νP)=∑NP=ην(μP-NPνP-NP)αP(νP-NP)(1-αP)(μP-νP)·Wn-1S(NP).The selection phase of pathogenP is described by (see (4) from Section, now with index P)
(15)WnS(NP)=∑ν=NPNmax(νPNP)βPNP(1-βP)(νP-NP)·WnM(νP).The multiplication phase of immune effectorE (producing antibody a) and memory M is described by (see (9) and (10) from Section 2.4, now with indices E and M)
(16)WnM(νE,νM)=∑NE=0Nmax∑NM=0Nmax-NE∑kE=0μE-NE∑kEM=0μE-kE-NEWn-1S(NE,NM)·(μE-NEkE)αE(μE-kE-NE)(1-αE)kE·(μE-kE-NEkEM)εEMkEM(1-εEM)(μE-kE-kEM-NE)·(μM-NMkM)αM(μM-kM-NM)(1-αM)kM·(μM-kM-NMkME)εMEkME(1-εME)(μM-kM-kME-NM),
where the two remaining indices are given by
(17)kM=μE+μM-νE-νM-kE,kME=νE+kEM+kE-μE.
The cut-off conditions are
(18)μE=ρE·NEμM=ρM·NMifρE·NE+ρM·NM≤Nmax,μE=νadjμM=Nmax-νadjifρE·NE+ρM·NM>Nmax,μP=ρP·NPifρP·NP≤Nmax,μP=NmaxifρP·NP>Nmax.The selection phase of immune effectorE and memory M is described by (see (12) from Section 2.4, now with indices E and M)
(19)WnS(NE,NM)=∑νE=NENmax∑νM=NMNmax-νEWnM(νE,νM)·(νENE)βENE(1-βE)(νE-NE)·(νMNM)βMNM(1-βM)(νM-NM).Each individual has a probabilityβ, that is, βP(E-), βE(P-), and βM(P-), being selected to survive, it replicates by a copy-probability α=0.1, multiplication-factor ρ, that is, ρP(E-), ρE(P-), and ρM(P-) with a mutation probability ε, that is, εEM(P-) and εME(P-). Thus, the multiplication factor ρ, the error probability ε, and the surviving probability β are parameters in function of averages P-=∑NP=0NP,maxW(NP)·NP and E-=∑NE,NM=0NmaxW(NE,NM)·NE (omitting the indices S/M and n for simplicity) that form the coupling (step functions with threshold values, see captions of Figures 8, 2 and 4).In Figure9 the case is shown, where the pathogen is eliminated (probability of extinction W(0)=1.0 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 W(0) persists below 1.0).Figure 9
Computer result of a pathogenP (with antigen A) coupled through averages to an immune system consisting of an effector E (producing antigen receptor a) and memory M. Case where pathogen is eliminated. Glimpses at generation (a) n=13, (b) n=84 and (c) n=200. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for N): probability WnS(NP) of finding 0≤NP≤Nmax individuals of the pathogen P (green). Lower left (triangle graph): probability WnS(NE,NM) of finding NE individuals of the immune effector E (red) and NM individuals of the immune memory M (blue), (0≤NE+NM≤Nmax). Upper right: extinction probability W(0) as a function of generations n of pathogen (green), immune effector, and immune memory (violet). Lower right: average of pathogen P- (green), of immune effector E- (red), and of immune memory M- (blue) as a function of generations n. Parameter values: maximal total number Nmax=32; αP=0.1, αE=0.1, and αM=0.1. For average E-<TE=1.5 (below threshold value, immune effector E inactive): ρP=6, βP=0.96. For average E->TE=1.5 (above threshold value, immune effector E active): ρP=2, βP=0.65. For average P->TP=0.5 (below threshold value, pathogen P not seen by immune system): ρE=2, ρM=2, εEM=0.50, εME=0.01, βE=0.90, and βM=0.93. For average P->TP=0.5 (above threshold value, pathogen P seen by immune system): ρE=6, ρM=2, εEM=0.50, εME=0.50, βE=0.95, and βM=0.93.Figure 10
Computer result of a pathogenP (with antigen A) coupled through averages to an immune system consisting of an effector E (producing antigen receptor a) and memory M. Case where the pathogen is reappearing while the immune system response is low. Glimpse at generation n=100. Parameter values other than Figure 2: For average E->TE=1.5 (above threshold value, immune effector E active): βP=0.75. For average P-<TP=1.5 (below threshold value, pathogen P not seen by immune system): βE=0.85. For average  P->TP=1.5 (above threshold value, pathogen P seen by immune system): βE=0.90.


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

We consider a simple case (Figure11) of a pathogen with two alleles at an A-to-B genlocus (one epitope): the pathogen PA (or pathogen PB, resp.) expressing antigen A (or antigen B).        The multiplication phase of pathogens PA and PB is described by ((9), (10) from Section 2.4, now with indices PA and PB)
(20)WnM(νPA,νPB)=∑NPA=0Nmax∑NPB=0Nmax-NPA∑kPA=0μPA-NPA∑kPAPB=0μPA-kPA-NPAWn-1S(NPA,NPB).(μPA-NPAkPA)αPA(μPA-kPA-NPA)(1-αPA)kPA.(μPA-kPA-NPAkPAPB)εPAPBkPAPB(1-εPAPB)(μPA-kPA-kPAPB-NPA).(μPB-NPBkPB)αPB(μPB-kPB-NPB)(1-αPB)kPB.(μPB-kPB-NPBkPBPA)εPBPAkPBPA(1-εPBPA)(μPB-kPB-kPBPA-NPB),
where the two remaining indices are given by
(21)kPB=μPA+μPB-νPA-νPB-kPA,kPBPA=νPA+kPAPB+kPA-μPA.Figure 11
Immune response to a variable pathogen (pathogen strains each with different antigens which are presented to the immune system). Two alleles at aAB-genlocus of the pathogen expressing antigen A or B: the immune effectors respond specifically (lock and key principle) by proliferation and producing antigen-receptor (antibody) a or b.The selection phase of pathogensPA and PB is described by ((12) from section 2.4, now with indices PA and PB)
(22)WnS(NPA,NPB)=∑νPA=NPANmax∑νPB=NPBNmax-νPAWnM(νPA,νPB).(νPANPA)βPANPA(1-βPA)(νPA-NPA).(νPBNPB)βPBNPB(1-βPB)(νPB-NPB).The immune effectorEa (or immune effector Eb) responding specifically, thus producing antigen-receptor a (or antigen-receptor b) which recognizes the antigen A (or antigen B) and eliminates pathogen PA (or pathogen PB, resp.). While the antigen “A” is the “lock” and the antigen-receptor “a” is the corresponding “key” (or the antigen B is the lock and the antigen-receptor b is another, but corresponding key.In Figure12 we show a computer result of a varying pathogen (PA with antigen A changes into PB with antigen B with a certain probability εPAPB,εPBPA vice versa, see equations (20)–(22) with indices PA und PB describing pathogens A and B) coupled through average values to an immune system against antigen A (and against antigen B) consisting of an effector Ea and memory Ma (and an effector Eb and memory Mb). There is no “a to b” or “b to a”—transition within the immune system. The pathogen expressing antigen A is nearly eradicated, but the mutant pathogen strain-expressing antigen B 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.Figure 12
Computer result of a varying pathogen (PA with antigen A changing into PB with antigen B) coupled through averages to an immune system against antigen A consisting of an effector Ea and memory Ma and against antigen B consisting of an effector Eb and memory Mb. Partial change of pathogen PA to pathogen PB escaping immune effector Ea and with delayed immune response of effector Eb. Glimpses at generation (a) n=22, (b) n=100. Upper left (triangle graph for 0 and scales 0.15 for W and 32 for N): probability WnS(NPA,NPB) of finding NPA individuals of the pathogen PA with antigen A and of finding NPB individuals of the pathogen PB with antigen B (green) (0≤NPA+NPB≤Nmax). Middle left (triangle graph): probability WnS(NEa,NMa) of finding NEa individuals of the immune effector (red) and NMa individuals of immune memory (blue) against antigen A (0≤NEa+NMa≤Nmax). Lower left (triangle graph): probability WnS(NEb,NMb) of finding NEb individuals of the immune effector (red) and NMb individuals of immune memory (blue) against antigen B (0≤NEb+NMb≤Nmax). Upper right: extinction probability W(0) as a function of generations n of pathogen (green), immune effector & memory against antigen A and immune effector & memory against antigen B (both violet). Middle and lower right: average of pathogen P- (green), of immune effector E- (red) and of immune memory M- (blue) against antigen A (middle right) and against antigen B (lower right) as a function of generations n. Parameter values others than Figure 2: average TPA=1.5  and average TPB=1.5, respectively (threshold value to switch immune system); εPAPB=0.01, εPBPA=0.001 (pathogen variability).


## 4. Maturation of T-Lymphocytes

As mentioned in Section1.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 LNaive, lymphocyte LSelf with strong affinity to self-peptides, matured lymphocyte LMat->I with weak affinity to foreign-peptides, this lymphocyte gets inactivated, matured lymphocyte LMat->A with strong affinity to foreign-peptides, this lymphocyte gets activated (Mat = matured, I = inactivated, A = activated). The multiplication phase is described by
(23)WnM(νLNaive,νLSelf,νLMat->I,νLMat->A)2=∑NLNaive=0Nmax∑NLSelf=0Nmax-NLNaive∑NLMat->I=0Nmax-NLNaive-NLSelf∑NLMat->A=0Nmax-NLNaive-NLSelf-NLMat->I22Wn-1S(NLNaive,NLSelf,NLMat->I,NLMat->A)22·(μLNaive-NLNaivekLNaive)αLNaive(μLNaive-kLNaive-NLNaive)(1-αLNaive)kLNaive22·(μLSelf-NLSelfkLSelf)αLSelf(μLSelf-kLSelf-NLSelf)(1-αLSelf)kLSelf22·(μLMat->I-NLMat->IkLMat->I)αLMat->I(μLMat->I-kLMat->I-NLMat->I)22·(1-αLMat->I)kLMat->I22·(μLMat->A-NLMat->AkLMat->A)αLMat->A(μLMat->A-kLMat->A-NLMat->A)22·(1-αLMat->A)kLMat->A,
where the remaining indices are given by
(24)kLNaive=μLNaive-νLNaive,kLSelf=μLSelf-νLSelf,kLMat->I=μLMat->I-νLMat->I,kLMat->A=μLMat->A-νLMat->A.
For the phase (a) in Figure 13, the multiplication phase is described by inserting formula (25) into formula (23)
(25)∑kLSelf=0μLSelf-NLSelf∑kLMat->I=0μLMat->I-NLMat->I16∑Permutation(LSelf,LMat->I,LMat->A)[(μLNaive-NLNaivekLNaiveLSelf)2·εLNaiveLSelfkLNaiveLSelf(1-εLNaiveLSelf)(μLNaive-NLNaive-kLNaiveLSelf)2·(μLNaive-NLNaive-kLNaiveLSelfkLNaiveLMat->I)2·εLNaiveLMat->IkLNaiveLMat->I2·(1-εLNaiveLMat->I)(μLNaive-NLNaive-kLNaiveLSelf-kLNaiveLMat->I)2·(μLNaive-NLNaive-kLNaiveLSelf-kLNaiveLMat->IkLNaiveLMat->A)2·εLNaiveLMat→AkLNaiveLMat→A2·(1-εLNaiveLMat->A)(μLNaive-NLNaive-kLNaiveLSelf-kLNaiveLMat->I-kLNaiveLMat->A)(μLNaive-NLNaivekLNaiveLSelf)],
where the remaining indices instead of (24) are given by
(26)kLNaiveLSelf=-μLSelf-νLSelf+kLSelf,kLNaiveLMat->I=-μLMat->I-νLMat->I+kLMat->I,kLNaiveLMat->A=μLNaive+μLSelf+μLMat->I-νLNaive-νLSelf-νLMat->I-kLSelf-kLMat->I,kLMat->A=μLNaive+μLSelf+μLMat->I+μLMat->A-νLNaive-νLSelf-νLMat->I-νLMat->A-kLSelf-kLMat->I.Figure 13
Computer result of the maturation process of T-lymphocytes. First a naïve T-lymphocyte (LNaive, green) in bone marrow or thymus undergoes T-cell receptor rearrangement (β-selection). T-cells with high affinity to self-peptides MHC (LSelf, black) 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 (LMat->I, blue and LMat->A, red) 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 (LMat->A, red), transforms into an effector cell and proliferates.  Glimpses at generation at  n=27 and n=83. Left: probabilities WnS(N) of finding N individuals of T lymphocytes (0≤N≤Nmax). Upper right: extinction probability W(0) as a function of generations n. Lower right: average of  T lymphocytes L- as a function of generations n. Parameter values and their change during the dynamics (a) n<50 (b) 50≤n<100 (c) 100≤n: maximal total number Nmax=16; αLNaive=0.1, αLSelf=0.1, αLMat->I=0.1, αLMat->A=0.1, ρLNaive=2, ρLSelf=1, ρLMat->I=1, ρLMat->A=1/4/1, βLNaive=0.95, βLSelf=0.75, βLMat->I=0.99/0.75/0.75, βLMat->A=0.99/0.98/0.75, εLNaiveLSelf=0.7, εLNaiveLMat->I=0.7, εLNaiveLMat->A=0.7.The selection phase of  lymphocytes is described by(27)WnS(NLNaive,NLSelf,NLMat->I,NLMat->A)2=∑νLNaive=NLNaiveNmax∑νLSelf=NLSelfNmax-νLNaive∑νLMat->I=NLMat->INmax-νLNaive-νLSelf∑νLMat->A=NLMat->ANmax-νLNaive-νLSelf-νLMat->I222WnM(νLNaive,νLSelf,νLMat->I,νLMat->A)222·(νLNaiveNLNaive)βLNaiveNLNaive(1-βLNaive)(νLNaive-NLNaive)222·(νLSelfNLSelf)βLSelfNLSelf(1-βLSelf)(νLSelf-NLSelf)222·(νLMat->INLMat->I)βLMat->INLMat->I(1-βLMat->I)(νLMat->I-NLMat->I)222·(νLMat->ANLMat->A)βLMat->ANLMat->A(1-βLMat->A)(νLMat->A-NLMat->A).


## 5. Discussion and Conclusion

Understanding the dynamics of both an invading pathogen and the response of the host’s immune system is an essential task in one’s attempt to positively influence the immune response of the given host. However, one already experiences difficulties in modeling the behavior of a single biological cell. A cell (as an element of a population of such cells) divides more frequently within a favorable environment and may form new strains by occasional errors in the copying process. It also dies more probably within a less favorable environment.How should one model such a step-by-step Darwinian process? Some may opt for numerically solving a set of first order differential equations, where time is continuous, and then examine the mainly nonlinear properties of variables (which represent large population rescaling). In contrast, we presented here a simple stochastic model of an entity undergoing a continued step-by-step Darwinian process, which is subdivided into two phases of multiplication (with variation) and selection. We describe this stochastic mathematically by a recursion formula (Galton-Watson type) for each phase, the dynamics of the system being evaluated numerically by computer, where the number of generation is an integer time-variable. The form of probability distributionW(N) changes in this system dynamics (with N being the number of individuals, and including the probability W(0)  of extinction). This is a great advantage of this approach.In addition, at a more fundamental level, one can suggest the following experiment to verify the Galton-Watson type dynamical stochastic process without mutation, described by equations (3) and (4) resulting in Figure 4 (case Section 2.3), or with copying errors leading to one mutant, described by equations (9)–(12) resulting in Figure 7 (case Section 2.4) and its parameter range: one prepares a steady-state condition of a bacterium-culture on a growing medium, where an antibiotic is added to the nutrient solution in a sublethal concentration, by repeated consecutive single-cell inoculation procedures. Then one can count bacteria by stopwatch the final single-cell inoculations carried out in parallel with the same nutrient solution (case Section 2.3) or with the nutrient solution charged additionally with another antibiotic of sub-lethal concentration (case Section 2.4) and plot the resultant time-dependant histogram.In this paper, we studied three types of behavior by analyzing both the pathogen and the host’s immune reaction with the proposed model system: (i) lasting pathogen elimination with buildup of immune memory, (ii) an oscillatory chronic case, where the pathogen is almost eliminated by the activated immune system, while during the subsequent relaxation of the immune system the pathogen is flaring up again, and (iii) the two-strain case, where the pathogen can vary its antigen at one epitope resulting in a prolonged immune-response.In order to map such a simple mathematical model of the immune-response to a real system, for example, a specific host, a specific pathogen, and a specific pathway, further work should consider the particular properties (e.g., the relative doubling rate) of the pathogen and the particular properties of the T  and B lymphocytes and other host properties as done by the aforementioned authors [2–4]. Explicitly considering genotype and phenotype should also be fruitful. Finally, one can find possible applications of the model in HIV, LCMV, influenza virus, herpes virus, mycobacterium tuberculosis, and plasmodium or trypanosomes. Supplementary material provided on the Website of CMMM available online at doi:10.1155/2012/784512.



---

*Source: 784512-2012-12-02.xml*