Dynamic Programming-Based Approach to Model Antigen-Driven Immune Repertoire Synthesis

Abstract

This paper presents a novel approach to modeling the repertoire of the immune system and its adaptation in response to the evolutionary dynamics of pathogens associated with their genetic variability. It is based on application of a dynamic programming-based framework to model the antigen-driven immune repertoire synthesis. The processes of formation of new receptor specificity of lymphocytes (the growth of their affinity during maturation) are described by an ordinary differential equation (ODE) with a piecewise-constant right-hand side. Optimal control synthesis is based on the solution of the Hamilton–Jacobi–Bellman equation implementing the dynamic programming approach for controlling Gaussian random processes generated by a stochastic differential equation (SDE) with the noise in the form of the Wiener process. The proposed description of the clonal repertoire of the immune system allows us to introduce an integral characteristic of the immune repertoire completeness or the integrative fitness of the whole immune system. The quantitative index for characterizing the immune system fitness is analytically derived using the Feynman–Kac–Kolmogorov equation.

1. Introduction

The main function of the human and animal immune system is to control and maintain antigenic homeostasis of the internal environment of the body. The realization of this function is associated with the activity of innate immune cells reacting to conservative antigenic structures of pathogens, and the reactions of cells of the adaptive immune system reacting with specific antigens of viruses and bacteria. Adaptive immune cells, B- and T-lymphocytes, have receptors on their surface that recognize a fragment of antigens in native form and in the form of complexes with receptors encoded by Major Histocompatibility Complex (MHC) genes (the main histocompatibility complex), respectively [1]. The effectiveness of the adaptive immune system is determined by the presence and number of clones of B- and T-lymphocytes expressing receptors that recognize antigens, i.e., B-cell receptors (BCR) and T-cell receptors (TCR), respectively. The formation of the clonal repertoire and its adaptive dynamics are realized via a sequence of division and differentiation stages in primary and secondary lymphoid organs. Most of the existing mathematical models describing the immune response development consider the reaction of a single lymphocyte clone with the most complementary receptor to the antigen (monoclonal dynamics) in accordance with the clonal selection principle.

The sets of antigens expressed by pathogen populations interacting at any given time with the human immune system are random samples from some general antigen population represented by 1407 pathogens capable of infecting humans [2]. Each clone of the lymphocytes expresses receptors specific to a specific antigen, while the clonal repertoire can also be considered as a limited sample of more than ${10}^{20}$ receptors of the general clonal population [3]. The totality of all clones makes up the repertoire of the immune system. The human immune system consists of about ${10}^{12}$ lymphocytes, forming up to ${10}^8$ clones, differing in the specificity of the antigen-recognizing receptors expressed by them. It is known that the distribution of clone numbers $C_i,i=1,\dots ,{10}^8$, follows a power law $p\left(C\right)\propto C^{-\alpha },\alpha \ge {\alpha }_{min}$ [4].

A variety of antigen-recognizing receptor structures of lymphocytes are formed during the development of lymphocytes, creating the basis for the selection of the functional repertoire of mature lymphocytes. The formation of a variety of antigenic receptors occurs as a result of a somatic recombination of innate gene segments encoding variable receptor regions (see for details [1]). A random rearrangement of the genes of developing lymphocytes leads to the generation of antigenic receptors of various specificities, passing further through the filters of positive and negative selection. The binding strength of the antigen-binding site of the receptor and the antigen epitope is called affinity and is expressed by the value of the dissociation constant in molar units. The main repertoire of antibodies has characteristic constants $K_d={10}^{-6}÷{10}^{-9}M$ and can improve by two orders of magnitude during antigen-induced maturation under the action of mutation, recombination, and selection processes in primary and secondary lymphoid organs.

Consider the sets (categories) of antigens involved in the formation of the receptor repertoire. At the initial selection stage taking place in the primary lymphoid organs, the sets of the spectrum of auto-antigens differ significantly from the foreign antigens shaping immune reactions in the secondary lymphoid organs. Let the number of elements of such sets be countable and equal to $i=1,2,\dots ,n$. Each antigen can be characterized (i.e., subject to coordinatization according to Weyl [5]) by some quantitative variable $r_i$, which may relate to the stereochemical structure of the antigen, coding genes, the probability of existence in the set under consideration, or the growth/fitness index of the antigen (let us call it a phenotype). Every pathogen ($i,i=1,2,\dots ,n$) expresses a number of antigens ${\mathbf{r}}_i$ which are recognized by the antigen-specific antibodies and receptors ${\mathbf{R}}_j,j=1,2,\dots ,m$ expressed on the immune cells. As noted before, the match (functional efficacy) between the immune receptor and a particular antigen is biochemically characterized by its affinity, i.e., the ratio of rate constants of association and dissociation. The perfect match corresponds to a $Key and Lock$ view of the antigen-receptor interaction.

The modeling of immune cell repertoire evolution remains an underexplored area. To describe the antibody diversity, maximum entropy models of sequence repertoire were proposed in [6]. The models predicted that the repertoire decomposes in several clusters and the sequence distribution obeys Zipf’s law. Importantly, the results suggested that the diversity of antibodies was dependent on both the genome-encoded sequences and the antigen-driven adaptation. Later on, a mathematical model of neutrally evolving pathogen and adaptation of the B-cell repertoire was described by a discrete random walk equation while the clonal immune responses were described by a discrete model of cell proliferation and death depending on the extent of cross-reactivity [7]. Assuming that the immune system follows an optimal strategy based on maximizing the antigenic space coverage at a minimal short-term metabolic cost, the model predicted a power-law distribution of the clonal sizes and specified conditions, in which the de novo synthesis of immune clone is more efficient in controlling mutating pathogen rather than relying on highly specific memory cells. The repertoire dynamics of T-cells in humans has been recently examined in [8]. To this end, a mathematical model of stochastic clonal dynamics by a geometric Brownian motion was formulated. The model was used to assimilate longitudinal repertoire sequencing data with characteristic time-scale of T-cell repertoire dynamics. For a comprehensive and critical review of mathematical models formulated to examine the population dynamics of immune repertoire, we refer to [9]. It is noted that the evolution of individual clones is described by Markovian dynamics with the expansion and death rates being dependent on some function which parameterizes the effects of specific antigens’ abundance and the competition between clones. An approach to model the dynamics of T cell clonotypes in response to virus infections for influenza and coronaviruses was proposed in [10]. The dynamics of T-cells’ clones is affected by a number of events (such as division, death, etc.) represented as a stochastic Markov process. A special focus is on cross-reactivity and competition between T-cells’ clones under homeostatic conditions and in response to the above-mentioned virus infections. Integrated analysis of immune repertoires is a fundamental issue in the analysis of amino acid sequencing-produced large datasets of T-cell receptors, Major Histocompatibilty Complexes, antigens, and antibodies. Unifying quantitative metrics of repertoire diversity and amino acid patterning have been proposed in [11]. Two powerful concepts from information theory, namely Shannon entropy and mutual information, are used to characterize the diversity in a population under study. The mean-field model is proposed in [12] to describe the homeostatic dynamics of naive T-cell clones under the impact of birth–death–thymic immigration processes. It was predicted that the clonal repertoire diversity is less dependent on heterogeneity in immigration rates but turns out to be affected by small variations in proliferation rates.

The general view of the corresponding models of the joint evolution of the viral population and B-lymphocytes can be represented by some generalized equations (see [13]) for the density distributions of antigens in the space of antigenic determinants $\mathbf{r}$, $V(t,\mathbf{r})$, and B-lymphocytes in the space antigen receptors $\mathbf{R}$, $B(t,\mathbf{R})$. In this study, we consider a fundamentally different approach to modeling the evolution of the clonal repertoire of the immune system. It is based on the ensemble dynamics of random trajectories of the evolutionary process represented by the Feynman-Kac–Kolmogorov equations and the Hamilton–Jacobi–Bellman dynamic programming methods [14], developed earlier in [15] and applied recently in a different context [16].

The models developed in this study focus on qualitative and numerical investigation into the joint dynamics patterns of pathogens and antigen-specific lymphocyte receptors in a certain metric space of the stereochemical structure (so-called ‘shape-space’) of expressed antigens (epitopes) and lymphocyte receptors in terms of the power of corresponding sets: ${\mathbf{r}}_i^{*}\left(t\right)$$i=1,2,\dots ,n$, and ${\mathbf{R}}_j^{*}\left(t\right)$, $j=1,2,\dots m$, $t\in [0,T]$, i.e., the numbers of actual epitopes and repertoire diversity. The processes of mutation and recombination can lead to a shift in the characteristics of the structure within certain spheres ${\mathcal{B}}_{\Delta }\left({\mathbf{r}}_i^{*}\right)$ and ${\mathcal{B}}_{\Delta }\left({\mathbf{R}}_j^{*}\right)$, respectively. For subsequent analysis, we will simplify the description of the position of epitopes and complementary receptors in the shape space by projections on one-dimensional axes for antigens and receptors, respectively.

In Section 2, we formulate the stochastic model of antigen-driven complementary lymphocyte receptor formation using the dynamic programming approach. In Section 3, the notions of fitness for antigen and immune system are defined. In Section 4, the immune-mediated control of replicating pathogens is examined. We conclude the analysis in Section 5.

2. Model of Antigen-Driven Lymphocyte Receptor Synthesis

2.1. Basic Model of Random Antigen Dynamics

We represent antigens in the form of Gaussian random processes $r_i\left(t\right)$, $i=1,2,\dots ,n$, $t\in [0,T]$ with the mathematical expectation ${\rho }_i$ and the correlation function, which is given by the formula

$$K_{r_i}=\mathbb{E}[(r_i\left(t\right)-{\rho }_i)\left(r_i(t+\gamma )\right)-{\rho }_i)]={\sigma }_i^2{e}^{-k_i\left|\gamma \right|},$$

(1)

where $k_i\ge 0$ is the decay parameter [17]. Each such process is generated by a stochastic differential equation of the following form [18]:

$$d{r}_i\left(t\right)=-k_i(r_i\left(t\right)-{\rho }_i)dt+{\sigma }_i\sqrt{2k_i}d{w}_t^i,i=1,2,\dots ,n.$$

(2)

Here, $w_t^i$ are independent Wiener random processes, and $k_i(r_i\left(t\right)-{\rho }_i)$ is the value of the drift of a random process from the reference value. The parameters $k_i,{\rho }_i,{\sigma }_i$ represent characteristics of randomly changing properties of antigens, realized in time as trajectories of Gaussian signals with noise in the form of a Wiener process. At the same time, ${\rho }_i$ characterizes the reference/stable value of the antigenic trait, ${\sigma }_i$—the intensity of the process of random variability (mutations and recombination), $k_i$—the characteristic rate of recovery (backward mutations) to the reference value, which can be associated with the effect of selection due to a higher fitness of the reference antigen type. In Figure 1, an example of random evolutionary dynamics for ten different antigens is presented.

2.2. Model of Antigen-Driven Lymphocyte Receptor Formation

Here, we introduce a description of the development of clones of lymphocytes with receptors complementary to the corresponding antigens, $R_i\left(t\right)$. To this end, consider the set of functions $\mathbf{r}\left(t\right)\doteq {[r_1\left(t\right),r_2\left(t\right),\dots ,r_n\left(t\right)]}^T$, $\mathbf{R}\left(t\right)\doteq {[R_1\left(t\right),R_2\left(t\right),\dots ,R_n\left(t\right)]}^T$ and controls $u_i\left(t\right)$, $i=1,2,\dots ,n$ describing the lymphocyte clone generation and selection processes. The proposed differential equations are proposed to take the following form:

$${\displaystyle \frac{d}{dt}}{R}_i\left(t\right)=u_i\left(t\right),R_i\left(0\right)=R_i^0,\left|u_i\left(t\right)\right|\le M_i,i=1,2,\dots ,n,M_i=\mathrm{const}>0.$$

(3)

Let us assume that the clonal repertoire of the immune system is successfully formed if there is such a set of functions $R_i\left(t\right)$ and controls $u_i\left(t\right)$, $i=1,2,\dots ,n$, for which the mathematical expectation of some norm of the mismatch function $W(\mathbf{r},\mathbf{R},t)$ between the receptors $R_i\left(t\right)$ and antigens $r_i\left(t\right)$ at some time, T, reaches a minimum value, e.g., the expectation of the squared $L_2$ norm:

$$W{(\mathbf{r},\mathbf{R},t)}_{t=T}={\mathbb{E}[\parallel \mathbf{r}\left(T\right)-\mathbf{R}\left(T\right)\parallel }_2^2].$$

(4)

2.3. Control Function for Synthesis of a Single Receptor

Consider the problem of finding control functions $u_i\left(t\right)$ and $R_i\left(t\right)$, $i=1,2,\dots ,n,$$0\le t\le T$ for the formation of an immune repertoire that implements the fulfillment of the condition (4). In order to minimize the mathematical expectation of the functional (4), we follow a dynamic programming approach using the Hamilton–Jacobi–Bellman equation (HJB). Let us first consider the case of one random Gaussian signal $r\left(t\right)$

$$dr\left(t\right)=(-k(r\left(t\right)-\rho )+R\left(t\right))dt+\sigma \sqrt{2k}d{w}_t,$$

(5)

and the single function $R\left(t\right)$ for receptor formation. By $W(r,R,\tau )$, $\tau =T-t$ the minimum value of the functional (4) is denoted, which can be achieved in the process of successful formation of a receptor recognizing the antigen in the problem (3) and (4) at time $\tau =T$. The function $W(r,R,\tau )$, $\tau =T-t$, is the solution of the following Cauchy problem for the HJB equation:

$${\displaystyle \frac{\partial W}{\partial \tau }}=-k(r-\rho ){\displaystyle \frac{\partial W}{\partial r}}+\underset{\left|u\right|\le M}{min}{\displaystyle \frac{\partial W}{\partial R}}u+{\sigma }^{2}k{\displaystyle \frac{{\partial }^{2}W}{\partial r^2}},$$

(6)

with the initial (terminal) condition

$$W(r,R,0)={(r-R)}^2.$$

(7)

The synthesis of control $u\left(t\right)$ for the immune response is determined by the following equalities:

$$u(r,R,\tau )=\left\{\begin{array}{cc}-M,\hfill & W_R(r,R,\tau )>0\hfill \\ +M,\hfill & W_R(r,R,\tau )<0\hfill \\ \mathrm{not}\mathrm{defined},\hfill & W_R(r,R,\tau )=0.\hfill \end{array}\right.$$

Here, the subscripts of the W function mean the corresponding partial derivatives. Let $W_R>0$. In this case, the solution of the HJB problem (6) and (7) denoted by $W^{+}(r,R,\tau )$ is

$$\begin{array}{cc}\hfill W^{+}(r,R,\tau )& ={(Z(r,\tau )-(R-M\tau ))}^2+{\sigma }^2(1-e^{-2k\tau }),\hfill \\ \hfill Z(r,\tau )& =e^{-k\tau }(r-\rho )+\rho .\hfill \end{array}$$

(8)

Due to the condition $W_R>0$, the solution represented by (8) is valid only in the domain

$$D^{+}=\{r,R,\tau :R\ge Z(r,\tau )+M\tau \}.$$

(9)

Proceeding similarly in the case when $W_R<0$, we obtain the solution $W^{-}(r,R,\tau )$:

$$\begin{array}{cc}\hfill W^{-}(r,R,\tau )& ={(Z(r,\tau )-(R+M\tau ))}^2+{\sigma }^2(1-e^{-2k\tau }),\hfill \\ \hfill Z(r,\tau )& =e^{-k\tau }(r-\rho )+\rho ,\mathrm{in}\mathrm{domain}\hfill \end{array}$$

(10)

$$D^{-}=\{r,R,\tau :R\le Z(r,\tau )-M\tau \}.$$

(11)

Consider the objective function defined as follows:

$$W(r,R,\tau )=\left\{\begin{array}{cc}{W}^{+}(r,R,\tau ),\hfill & (r,R,\tau )\in D^{+}\hfill \\ W^{-}(r,R,\tau ),\hfill & (r,R,\tau )\in D^{-}\hfill \\ W^0(r,R,\tau ),\hfill & (r,R,\tau )\in D^0.\hfill \end{array}\right.$$

(12)

Here,

$$\begin{array}{cc}\hfill W^0(r,R,\tau )& =Z^2(r,\tau )+{\sigma }^2(1-e^{-2k\tau }),\hfill \\ \hfill Z(r,\tau ;\rho )& =e^{-k\tau }(r-\rho )+\rho ,\mathrm{in}\mathrm{domain}\hfill \end{array}$$

(13)

$$D^0=\{r,R,\tau :Z(r,\tau )-M\tau <R<Z(r,\tau )+M\tau \}.$$

(14)

Function $W(r,R,\tau )$ is the solution of the HJB Equation (6). This function is continuous with respect to the set of variables $r,R,\tau$, has continuous first and second derivatives with respect to the variable r, and continuous derivatives with respect to the variables R and $\tau$, and also satisfies conditions (12) in the domain $D^{+}$ and $D^{-}$. Hence, the synthesis of the control function for the formation of the clone with the receptor R can be defined by the following expressions, which depend on subsets determined adaptively by $r,R,\tau$:

$$u(r,R,\tau ;\rho )=\left\{\begin{array}{cc}-M,\hfill & R\ge Z(r,\tau )+M\tau ,D^{+}\hfill \\ +M,\hfill & R\le Z(r,\tau )-M\tau ,D^{-}\hfill \\ 0,\hfill & Z(r,\tau )-M\tau <R<Z(r,\tau )+M\tau ,D^0.\hfill \end{array}\right.$$

(15)

The above three domains differing with respect to the mode of control synthesis of antigen-related receptor specificity are shown in Figure 2.

2.4. Control Function Accounting for the Gaussian Variability of Antigen Dynamics

The random process variable $r\left(t\right)$ at every time t follows a Gaussian law because the empirical values of the parameters $\rho$, $\sigma$ and k of the underlying governing Equation (5) represent random realizations of some Gaussian distributions. Therefore, it is necessary to perform an averaging of the mismatch function ($W(r,R,\tau )$) over all realization of the Gaussian distribution $r\left(t\right)$:

$$\begin{array}{ccc}\hfill \overline{W}(R,\tau )& =& {\int }_{-\infty }^{+\infty }W(r,R,\tau )P_{\sigma }(r-\rho )dr,\hfill \\ \hfill P_{\sigma }(r-\rho )& =& {\displaystyle \frac{1}{2\sqrt{\pi }\sigma }}{e}^{-{\displaystyle \frac{{(r-\rho )}^2}{2{\sigma }^2}}}.\hfill \end{array}$$

(16)

Following analysis similar to the previous subsection, one can show that, for the following expressions, the averaged function can be derived:

$$\overline{W}(R,\tau )=\left\{\begin{array}{cc}{(\rho -(R-M\tau ))}^2+{\sigma }^2,\hfill & R>\rho +M\tau \hfill \\ {(\rho -(R+M\tau ))}^2+{\sigma }^2,\hfill & R<\rho -M\tau \hfill \\ {(\rho -R)}^2+{\sigma }^2,\hfill & \rho -M\tau \le R\le \rho +M\tau .\hfill \end{array}\right.$$

(17)

For the above defined averaged functional $\overline{W}(R,\tau )$, the HJB Equation (6) takes the form of a hyperbolic partial differential equation (PDE):

$${\displaystyle \frac{\partial \overline{W}(R,\tau )}{\partial \tau }}=\underset{\left|u\right|\le M}{min}\left\{{\displaystyle \frac{\partial \overline{W}}{\partial R}}(R,\tau )u\right\},$$

(18)

with the initial (terminal) condition

$$\overline{W}(R,0)={\sigma }^2+{(R-\rho )}^2,$$

(19)

in the domains ${\overline{D}}^{+}=\{{\displaystyle \frac{\partial \overline{W}}{\partial R}}\ge 0\}=\{R,\tau :R\ge \rho +M\tau \}$, and ${\overline{D}}^{-}=\{{\displaystyle \frac{\partial \overline{W}}{\partial R}}\le 0\}=\{R,\tau :R\le \rho -M\tau \}$.

In domain ${\overline{D}}^0=\{R,\tau :\rho -M\tau <R<\rho +M\tau \}$, the initial value problem for HJB equation reduces to

$${\displaystyle \frac{\partial \overline{W}(R,\tau )}{\partial \tau }}=0,\overline{W}(R,0)={\sigma }^2+{(R-\rho )}^2.$$

(20)

Hence, $\overline{W}(R,\tau )={\sigma }^2+{(R-\rho )}^2$ in ${\overline{D}}^0$, if the $\tau$-axis belongs to this domain.

Synthesis of control function for the immune repertoire formation in relation to R, $\tau$ and $\rho$ can be derived by specifying the derivatives of the functions in (17) by variable R. It takes the following form:

$$\overline{U}(R,\tau )=\left\{\begin{array}{cc}-M,\hfill & \phantom{\rho -M\tau \le }R\ge \rho +M\tau ,{\overline{D}}^{+}\hfill \\ +M,\hfill & \phantom{\rho -M\tau \le }R\le \rho -M\tau ,{\overline{D}}^{-}\hfill \\ 0,\hfill & \rho -M\tau <R<\rho +M\tau ,{\overline{D}}^0.\hfill \end{array}\right.$$

(21)

Consider the evolution of the function $R\left(t\right)$ under the action of controls (21) with the initial condition $R(t=0)=R^0$ (note that this time value t corresponds to the initial value $\tau =T$ in backward time):

$${\displaystyle \frac{d}{dt}}R\left(t\right)=\overline{U}(R,\tau ;\rho ).$$

(22)

If $R>\rho +M\tau$, then $\overline{U}(R,\tau ;\rho )=-M$, and the solution of the problem (22) is $R\left(\tau \right)=-M\tau +R^0$. By virtue of the last condition, $R=-Mt+R^0>\rho +M\tau$. Therefore, the condition $R^0>\rho +MT$ needs to be satisfied.

Following a similar consideration when $R<\rho -M\tau$, and $\overline{U}(R,\tau ;\rho )=+M$, we obtain that the condition to be satisfied is $R^0<\rho -MT$. Overall, depending on the initial condition $R^0$, the function $R\left(\tau \right)$ is determined as follows:

$$R\left(t\right)=\left\{\begin{array}{cc}-Mt+R^0,\hfill & \phantom{\rho -MT\le }{R}^0\ge \rho +MT\hfill \\ +Mt+R^0,\hfill & \phantom{\rho -MT\le }{R}^0\le \rho -MT\hfill \\ R^0,\hfill & \rho -MT<R^0<\rho +MT,\hfill \end{array}\right.$$

(23)

whereas the cost functional $\overline{W}(R\left(t\right),T)$ at time T takes the following values:

$$\overline{W}(R\left(t\right),T)=\left\{\begin{array}{cc}{(\rho -(R^0-MT))}^2+{\sigma }^2,\hfill & \phantom{\rho -MT\le }{R}^0\ge \rho +MT\hfill \\ {(\rho -(R^0+MT))}^2+{\sigma }^2,\hfill & \phantom{\rho -MT\le }{R}^0\le \rho -MT\hfill \\ {(\rho -R^0)}^2+{\sigma }^2,\hfill & \rho -MT<R^0<\rho +MT.\hfill \end{array}\right.$$

(24)

3. Fitness Characterization of the Adapting Immune System

3.1. Complementarity to Antigen

The difference $\delta \overline{W}=\overline{W}(0,T)-\overline{W}(R\left(t\right),T)$ represents the amount by which the affinity of the (clonal) receptor to the antigen increases, and accordingly, the fitness of the antigen decreases. It can be quantified according to the following formulas:

$$\overline{\delta W}=\left\{\begin{array}{cc}(R^0-MT)(2\rho -(R^0-MT)),\hfill & \phantom{\rho -MT\le }{R}^0\ge \rho +MT\hfill \\ (R^0+MT)(2\rho -(R^0+MT)),\hfill & \phantom{\rho -MT\le }{R}^0\le \rho -MT\hfill \\ R^0(2\rho -R^0),\hfill & \rho -MT<R^0<\rho +MT.\hfill \end{array}\right.$$

(25)

In general, the antigen-receptor mismatch function $W(\mathbf{r},\mathbf{R},\tau )$, with $\mathbf{r}=(r_1,r_2,\dots ,r_n)$, $\mathbf{R}=(R_1,R_2,\dots ,R_n)$, $\tau =T-t$, can be represented as a sum of functions

$$W(\mathbf{r},\mathbf{R},\tau )=\sum _{i=1}^n{W}_i(r_i,R_i,\tau ),$$

(26)

where the functions $W_i(r_i,R_i,\tau )$ are the solutions of the backward in time Cauchy problem for the HJB equation:

$${\displaystyle \frac{\partial W_i}{\partial \tau }}=-k_i(r_i-{\rho }_i){\displaystyle \frac{\partial W_i}{\partial r_i}}+\underset{|u_i|\le M_i}{min}{\displaystyle \frac{\partial W_i}{\partial R_i}}{u}_i+{\sigma }_i^2{k}_i{\displaystyle \frac{{\partial }^2{W}_i}{\partial r_i^2}},$$

(27)

with the initial (terminal) condition

$$W_i(r_i,R_i,0)={(r_i-R_i)}^2,\tau =T-t.$$

(28)

One way to define the immune system fitness $F{t}_{IS}\left(t\right)$ (when the antigen replication is not considered) is the use of the function representing the squared mismatch between $\mathbf{r}$ and $\mathbf{R}$:

$$F{t}_{IS}\equiv \delta \overline{W}=\sum _{i=1}^n\delta \overline{W_i}=\sum _{i=1}^n(\overline{W_i}(0,T)-\overline{W_i}(R_i\left(t\right),T)).$$

(29)

Note that the synthesis of the adaptive repertoire in the immune system in the areas of $D_i^{+}=\{{\rho }_i,\tau :R_i\ge {\rho }_i+M_i\tau \}$, $D_i^{-}=\{{\rho }_i,\tau :R_i\le {\rho }_i-M_i\tau \}$, $D_i^0=\{{\rho }_i,\tau :{\rho }_i-M_i\tau <R_i<{\rho }_i+M_i\tau \}$, is described in a similar way to the previously presented.

3.2. Pathogen Load Control Based

Here, we develop a basic mathematical model for pathogen load ($v\left(t\right)$) dynamics with its antigenic trait parameter $r\left(t\right)$ under control of the adaptive immune response $R\left(t\right)$. We consider the following SDE describing the dynamics of antigenic parameter $r\left(t\right)$, which determines the antigenicity of pathogen $v\left(t\right)$:

$$dr\left(t\right)=-k(r\left(t\right)-\rho +R\left(t\right))dt+\sigma \sqrt{2k}d{w}_t,$$

(30)

Suppose that the antigenic properties of the pathogen set by the variable $r\left(t\right)$ are directly related to its fitness, i.e., affect the dynamics of the pathogen abundance $v\left(t\right)$ according to the Gompertz equation:

$$dv\left(t\right)=r\left(t\right)v\left(t\right)(S_v-log\left(v\left(t\right)\right))dt,S_v=\mathrm{const}>0.$$

(31)

Assuming $f\left(t\right)=log\left(v\right(t\left)\right)$, we obtain:

$$df\left(t\right)=r\left(t\right)(S-f(t\left)\right)dt,S=\mathrm{const}.$$

(32)

Consider a system of SDEs (31) and (32). The function $F(f,r,\tau ;R(\tau \left)\right),\tau =T-t$, representing the mathematical expectation of the size (log-transformed) of the pathogen population f at time $t=T$ is a solution to the Cauchy problem for the Feynman–Kac–Kolmogorov (FKK) equation:

$${\displaystyle \frac{\partial F}{\partial \tau }}=r(S-f){\displaystyle \frac{\partial F}{\partial f}}+[-k(r-\rho )+R\left(\tau \right)]{\displaystyle \frac{\partial F}{\partial r}}+{\sigma }^{2}k{\displaystyle \frac{{\partial }^{2}F}{\partial r^2}},$$

(33)

subject to the terminal condition

$$F(f,r,0;R(T\left)\right)=f;F(S,r,\tau ;R(\tau \left)\right)=S>0;\tau =T-t.$$

(34)

The solution to problems (33) and (34) can be presented in an explicit form as follows:

$$\begin{array}{cc}\hfill F(f,r,\tau ;R(\tau \left)\right)& =S-(S-f)exp(z_1(r,\tau )+z_2\left(\tau \right)),\hfill \\ \hfill z_1(r,\tau )& =-\rho \tau +{\displaystyle \frac{r-\rho }{k}}(e^{-k\tau }-1)+{\int }_0^{\tau }R\left(\tau \right)(1-e^{-k\tau })/kd\tau ,\hfill \\ \hfill z_2\left(\tau \right)& ={\displaystyle \frac{{\sigma }^2}{k}}{\int }_0^{\tau }{(e^{-k\tau }-1)}^{2}d\tau .\hfill \end{array}$$

(35)

This solution corresponds to a single process realization of the control function $u(r,R,\tau )$ provided by formula (15).

To take into account the Gaussian variability of the antigenic variable $r\left(t\right)$ at any given time t, we specify the solution to problem (33) and (34) using an averaged control function $\overline{U}(R,\tau )$ similar to expressions (21). To this end, it is necessary to calculate the mathematical expectation of the function $\overline{F}(f,r,\tau ;R\left(\tau \right))$, which represents the stochastic dynamics of the pathogen growth characteristic, for all Gaussian distributions of the value r:

$$\begin{array}{ccc}\hfill \overline{F}(f,\tau ;R\left(\tau \right))& =& {\int }_{-\infty }^{+\infty }F(f,r,\tau ;R\left(\tau \right))P_{\sigma }(r-\rho )dr,\hfill \\ \hfill P_{\sigma }(r-\rho )& =& {\displaystyle \frac{1}{2\sqrt{k\pi }\sigma }}{e}^{-{\displaystyle \frac{{(r-\rho )}^2}{2{\sigma }^2}}}.\hfill \end{array}$$

(36)

Making some transformations, we obtain the following representation

$$\begin{array}{cc}\hfill \overline{F}(f,\tau ;R\left(\tau \right))& =S-(S-f)exp(z_3\left(\tau \right)+z_4\left(\tau \right)),\hfill \end{array}$$

(37)

where

$$\begin{array}{cc}\hfill z_3\left(\tau \right)& =-\rho \tau +{\displaystyle \frac{{\sigma }^2}{2k}}{(e^{-k\tau }-1)}^2+{\displaystyle \frac{{\sigma }^2}{k}}{\int }_0^{\tau }{(e^{-k\tau }-1)}^{2}d\tau ,\hfill \\ \hfill z_4\left(\tau \right)& ={\int }_0^{\tau }R\left(\tau \right)(1-e^{-k\tau })/kd\tau .\hfill \end{array}$$

(38)

Now we can define the fitness of the immune system with respect to its ability to control the single antigen load. Considering the value $\overline{F}(f,\tau ;0)$, let us denote the function $\overline{F}(f,\tau ;R\left(\tau \right))$ with $R\left(\tau \right)=0$. It follows from expressions (37) and (38) that

$$\begin{array}{c}\hfill {\displaystyle \frac{S-\overline{F}(f,\tau ;0)}{S-\overline{F}(f,\tau ;R\left(\tau \right))}}=e^{-z_4\left(\tau \right)}.\end{array}$$

(39)

Note that the function $z_4\left(\tau \right)>0$. Consider $\overline{F}(f,\tau ;R\left(\tau \right))=\overline{F}(f,\tau ;0)-\delta F,\delta F>0$. The value $\delta F=(S-\overline{F}(f,\tau ;0))(e^{-z_4\left(\tau \right)}-1)$ represents the decrease in the abundance of the pathogen as a result of successful adaptation of the immune system repertoire. It can be considered as a measure of the fitness of the immune system with respect to a single pathogen.

Based on the above analysis, a global fitness of the immune system with respect to its ability to control multiple sets of varying pathogens can be derived. In the case of multiple pathogens, $F(\mathbf{r},\mathbf{R},\tau )$, where $\mathbf{r}=(r_1,r_2,\dots ,r_n)$, $\mathbf{R}=(R_1,R_2,\dots ,R_n)$, $\tau =T-t$, is a sum of single pathogen-related functions

$$F(\mathbf{r},\mathbf{R},\tau )=\sum _{i=1}^n{F}_i(r_i,\tau ;R_i).$$

(40)

In turn, the functions $F_i(r_i,\tau ;R_i)$ are the solutions to the backward in time Cauchy problem for the FKK equation:

$${\displaystyle \frac{\partial F_i}{\partial \tau }}=r_i(S_i-f_i){\displaystyle \frac{\partial F_i}{\partial f_i}}+[-k_i(r_i-{\rho }_i)+R_i\left(\tau \right)]{\displaystyle \frac{\partial F_i}{\partial r_i}}+{\sigma }_i^2{k}_i{\displaystyle \frac{{\partial }^2{F}_i}{\partial r_i^2}},$$

(41)

subject to the terminal condition

$$F_i(f_i,r_i,0;R\left(T\right))=f_i.$$

(42)

4. Numerical Simulations of Immune System–Pathogen Co-Adaptation

4.1. Repertoire Synthesis for a Single Pathogen

Here, we examine the quantitative dynamics of pathogen load under the control of the adapting immune system. In the Gompertz equation studied above, the fitness of the pathogen was assumed to be directly proportional to $r\left(t\right)$. Obviously, the synthesis of lymphocyte clone with the receptor R complementary to the specific antigen r should result in the reduction in the pathogen fitness and the net growth rate of the pathogen. To parameterize this effect, we assume that the pathogen growth rate depends on the difference between $r-R$ according to the following function:

$$\mu (r,R)=1-(1-\gamma )exp\left(-{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{|r-R|}{{\sigma }_{\Delta }}}\right)}^m\right),\gamma \in (0,1),{\sigma }_{\Delta }>0,m>0,$$

(43)

so that the dynamics of the pathogen load $v\left(t\right)$ is described by equation:

$${\displaystyle \frac{d}{dt}}v\left(t\right)=\beta {\displaystyle \frac{r\left(t\right)}{\rho }}\mu \left(r\right(t),R(t\left)\right)v\left(t\right)(S_v-log\left(v\left(t\right)\right)),$$

(44)

where the antigen-specific function $r\left(t\right)$ is governed by the SDE equation:

$$dr\left(t\right)=-k(r\left(t\right)-\rho )dt+\sigma \sqrt{2k}d{w}_t.$$

(45)

The immune system adaptation, i.e., the complementary clone $R\left(t\right)$ formation, is described by the ODE with a piecewise-constant right-hand side:

$${\displaystyle \frac{d}{dt}}R\left(t\right)=\overline{U}(R,t;\rho ,M,T),t\in [0,T],$$

(46)

where the synthesis of control function $\overline{U}(R,t;\rho ,M,T)$ is defined in (21). The above parameterization $\mu (r,R)$ of the fitness dependence on the complementarity of the lymphocyte receptor suggests that the replication rate decreases from $\beta {\displaystyle \frac{r\left(t\right)}{\rho }}$ (no adaptation) to $\gamma ·\beta {\displaystyle \frac{r\left(t\right)}{\rho }}$ (fully adapted). Importantly, the three parameters of the fitness function, i.e., $\gamma$, ${\sigma }_{\Delta }$ and m, can be interpreted as the pathogen growth inhibition (no inhibition when $\gamma =1$ and full inhibition for $\gamma =0$), the cross-reactivity/specificity and the sensitivity threshold, respectively.

To obtain the numerical realizations of the model (44)–(46), we used the SRIW1 solver for SDEs, which is an adaptive strong order 1.5 and weak order 2.0 in the Itô sense, from the DifferentialEquations.jl package in Julia language. The event detection for switching between the values of the piecewise-constant function $\overline{U}(R,t;\rho ,M,T)$, i.e., the localization of the domains ${\overline{D}}^{+}$, ${\overline{D}}^{-}$ and ${\overline{D}}^0$, was implemented via callback functions.

Representative evolution of the mutating pathogen (six individual realizations), the dynamics of synthesis of the complementary receptors, and the pathogen fitness reduced by adapting immune system are shown in Figure 3.

An increase in the cross-reactivity/specificity of the synthesized clonal receptors associated with two-fold larger value of ${\sigma }_{\Delta }$ results in a slightly delayed growth of the pathogen according to the Gompertz equation. This is demonstrated in Figure 4. However, a two-fold decrease in the sensitivity threshold parameter m in the pathogen fitness function to the synthesized clonal receptors does not have a strong impact on the pathogen growth and the replication rate dynamics. This follows from comparing the trajectories in Figure 3 and Figure 5.

A strong impact on the pathogen dynamics results from five-fold increase in the pathogen growth inhibition associated with a five times decrease in the respective parameter $\gamma$. The respective control of the pathogen dynamics is shown in Figure 6. One can see that the fitness values are reduced much faster, and the pathogen loads are on average smaller. Notice that the effect depends on the initial growth rates randomly sampled from a Gaussian distribution, and it becomes stronger for the smaller values.

The rate of the synthesis of complementary receptor in the immune system is represented by parameter M in the respective control synthesis Equation (46). We examined the effect of a two-fold increase in M. The overall result was an earlier completion of the control process and hence, the start-time for efficient reduction in the pathogen fitness was shorter. This is demonstrated in Figure 7 as compared to Figure 3.

All pathogens differ in the rate of their antigenic evolution. In the presented model, it is associated with parameter $\sigma$. Its increase from the basic set value 0.3 to 1 results in a larger variation of the pathogen antigenicity, as shown in Figure 8. It shows that faster mutating antigens represent ‘moving targets’ in the antigenic space that are more difficult for the immune system to control. Indeed, one can see that the amplitude of fluctuations of the individual realizations of pathogen dynamics increase substantially, as well as the overall growth of the pathogen. An increase in the degree of pathogen inhibition in the situation of highly variable intensity of random antigen fluctuations results in a broader variation of the replication rate dynamics by increasing the lower values of the replication rates (see Figure 9). Consequently, some of the pathogen load dynamics are characterized by a much slower growth kinetics.

It is known that under certain conditions the antibodies and their producing B-cells with complementary antigen-recognizing receptors can lead to an enhancement rather than inhibition of associated infections [19,20,21]. The developed model allowed us to examine such modes of antigen-immune system interaction. To this end, we simulated the impact of a broad initial distribution of complementary receptors on pathogen dynamics. The results are presented in Figure 10. One can see that a substantial diversity in pathogen fitness takes place.

4.2. Repertoire Synthesis for Multiple Pathogens

The immune system faces a continuous forcing from multiple antigens. The ability to control them can be considered as an overall fitness of the immune system. Analytical derivation of the immune system fitness evolution via the synthesis of complementary lymphocyte receptor clones was examined in the previous section for a simple model with the antigen fitness related to the antigenicity value $r\left(t\right)$ only. However, the adaptation of the immune system changes the fitness of the pathogens, and we used the Gompertz model with parameterization of the growth rate as a function of the mismatch $\left|r\right(t)-R(t\left)\right|$ to computationally explore the effect of variation in antigen-recognizing clonal receptor parameters for one pathogen in the above subsection. Here, we trace simultaneous dynamics of multiple pathogens and the response of the immune system in terms of the synthesis of a spectrum of complementary receptors. To this end, the following set of equations was considered: (2), (41), (44) and (46). Figure 11 presents the characteristic features of the system dynamics. One can see that an increasing fitness of the immune system results in the inhibition of the dynamics of pathogen load (solid lines) and the relative abundance of specific pathogen load (% of the maximal) (Figure 11).

Finally, we examined the joint effect of the pathogen growth inhibition degree $\gamma$ and the cross-reactivity/sensitivity parameter ${\sigma }_{\Delta }$ on the degree of the pathogen load reduction due to optimal synthesis of complementary receptors. The results of tracing the mean cumulative number of pathogens as a function of $\gamma$ and ${\sigma }_{\Delta }$, i.e., the value of

$${\langle AU{C}_T\left(v\right)\rangle }_{\gamma ,{\sigma }_{\Delta }}={\displaystyle \frac{1}{N}}\sum _{i=1}^N\left({\int }_0^T{v}_i\left(t\right)dt\right),N=1000,T=30,$$

are presented in the heatmap-type Figure 12, expressed as a percentage of reduction from the uncontrolled case

$${\langle AU{C}_T\left(v\right)\rangle }_{\gamma =1,{\sigma }_{\Delta }=0}.$$

5. Conclusions

This paper presents a new approach to modeling the repertoire of the immune system and its adaptation in response to the evolutionary dynamics of pathogens associated with their genetic variability. The processes of the formation of new receptor specificity of lymphocytes (the growth of their affinity during maturation) are described by an ODE with piecewise-constant right-hand side. Optimal control synthesis is based on the solution of the Hamilton–Jacobi–Bellman equation implementing the dynamic programming approach for controlling the Gaussian random processes generated by SDE with the noise in the form of the Wiener process. The solution of the Cauchy problem exists and is analytically derived in our study. Whether it is a unique one is a special issue. It goes beyond the aim of our work and requires a systematic thorough analysis following the qualitative theory of PDEs [22].

The proposed description of the clonal repertoire of the immune system allows us to introduce an integral characteristic of the completeness of the repertoire or the integrative fitness of the whole system, i.e.,

$$F{t}_{IS}=1-{\displaystyle \frac{{\langle AU{C}_T\left(v\right)\rangle }_{\gamma ,{\sigma }_{\Delta }}}{{\langle AU{C}_T\left(v\right)\rangle }_{\gamma =1,{\sigma }_{\Delta }=0}}}.$$

Note that a quantitative index for characterizing the immune system fitness for a simple consideration of the pathogen fitness was analytically derived using the Feynman–Kac–Kolmogorov equation in Section 2.

We used a simple description of the pathogen dynamics in the form of the Gompertz equation; however, the immune system functions to eliminate the foreign antigens from the host organisms. A direct extension of the modeling approach presented in this study would be incorporation of the equations (see [13]) describing the activation and clonal expansion of the immune response $B(t,\mathbf{R})$ and immune-mediated elimination of the pathogen $V(t,\mathbf{r})$ (rather than a mere reduction in its intrinsic growth rate):

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial t}}V(t,\mathbf{r})& =\underset{Antigenicvariation}{\underbrace{\nabla ·\left(D_{rV}\nabla V(t,\mathbf{r})\right)}}\hfill \\ \hfill & +\underset{Infectiongrowth}{\underbrace{b_V\left(\mathbf{r}\right)V(t,\mathbf{r})\left(K_V-J\left(V\right)\right)}}-\underset{Naturaldeath}{\underbrace{d_{V}V(t,\mathbf{r})}}\hfill \\ \hfill & -\underset{Immune-mediatedelimination}{\underbrace{d_{IC}V(t,\mathbf{r}){\int }_{{\Omega }_R}I(\mathbf{r}-\mathbf{R})B(t,\mathbf{R})d\mathbf{R}}}+\underset{Randomforcing}{\underbrace{{\sigma }_V{\xi }_V(t,\mathbf{r}),}}\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial t}}B(t,\mathbf{R})& =\underset{Primaryreceptorselection}{\underbrace{S(t,\mathbf{R})}}+\underset{Receptorsomatichypermutation}{\underbrace{\nabla ·\left(D_{RB}\nabla B(t,\mathbf{R})\right)}}\hfill \\ \hfill & +\underset{Bell-shapedclonalgrowth}{\underbrace{b_B\left(\mathbf{R}\right)\varphi \left(V(t,•)\right)B(t,\mathbf{R})\left(K_B-J\left(B\right)\right)}}-\underset{Naturaldeath}{\underbrace{d_{B}B(t,\mathbf{R})}}\hfill \\ \hfill & -\underset{Activartion-inducedapoptosis}{\underbrace{d_{AB}C(t,\mathbf{R}){\int }_{{\Omega }_r}\psi \left(V(t,\mathbf{r})\right)A(\mathbf{r}-\mathbf{R})d\mathbf{r}}}+\underset{Randomforcing}{\underbrace{{\sigma }_B{\xi }_B(t,\mathbf{R})}}\hfill \\ \hfill & -\underset{Affinitymaturation}{\underbrace{{\nabla }_R·\left({\chi }_{B}C(t,\mathbf{R}){\nabla }_{r-R}V(t,\mathbf{r})\right)}}.\hfill \end{array}$$

(48)

Further research will be related to the application of the above novel class of mathematical models and the analyses techniques developed earlier in a different context (see [16]) for host immune system–pathogen co-evolution to study: (1) the evolutionary dynamics of immune processes in healthy hosts and during infectious diseases under the influence of antigenic forcing, taking into account the frequency and temporal characteristics of the appearance of pathogens; (2) the regulation of the immune system repertoire required for optimal implementation of the protection function; (3) implementation of selection/maturation processes of antigenic lymphocyte receptors for a proper parameterization of appropriate control synthesis; and (4) the treatment of inverse problems for determining the parameters of models of adaptive dynamics of the clonal repertoire on evolving fitness landscapes.