The role of T-helper/T-regulator lymphocyte ratio on the adaptive immune response: a mathematical model Alessia Annibale Department of Mathematics, King’s College London, The Strand, London WC2R 2LS, UK Institute for Mathematical and Molecular Biomedicine, King’s College London, Hodgkin 7 1 Building, London SE1 1UL, UK 0 2 n a Abstract formation from them. Another reason is that, J while experiments in vitro and in mice have pro- 6 In this work we use a combination of statistical vided tremendous insights into the immune sys- 2 physics and dynamicalsystems approaches, to an- tems, the extrapolation of results to humans is ] alyze the response to an antigen of a simplified not always obvious [30], so theoretical frameworks B modeloftheadaptiveimmunesystem,whichcom- are demandedto modeltheimmunesystem in hu- C prises B, T helper and T regulatory lymphocytes. mans, where experimental work is very limited. . o Results show that the model is remarkably robust i A major area of focus in immunology in re- b against changes in the kinetic parameters, noise cent years has been regulatory T cells, which sup- - q levels, and mechanisms that activate T regulatory press other immune cells. These are believed to [ lymphocytes. In contrast, the model is extremely have several important functions, including main- 1 sensitive to changes in the ratio between T helper tainance ofhomeostasis, immunological unrespon- v and T regulatory lymphocytes, exhibiting in par- 2 siveness to self-antigens, and prevention of au- ticular a phase transition, from a responsive to an 1 toimmune diseases [31], and in the last few years 7 immuno-suppressed phase, when the ratio is low- have witnessed an explosion of knowledge driven 7 ered below a critical value. To the best of our 0 by detailed experimental work [32, 33, 34]. knowledge, this is the first mathematical study to . 1 In addition, recent studies have suggested that support the validity of the T-helper/T-regulator 0 the imbalance between regulatory and conven- 7 lymphocyte ratio as an index of immunosuppres- 1 sion, with a potential for prognostic monitoring. tional T cells may be an index of immunosup- : pression in cancer patients [35, 36], and may v i play a crucial role in the progression of HIV in- X fections and response to antiretroviral therapies r a [37, 38, 39], immunosenescence [40] and in inflam- 1 Introduction mation and autoimmune diseases [41]. However, the lack of reliable markers for distinguishing reg- Recent years have seen a surge of mathemat- ulatory from conventioanl T cells, and variations ical and computational models of the immune in the total lymphocyte count make regulatory T system, ranging from ordinary differential equa- cells counting controversial [36]and ageneral con- tions approaches [1, 2, 3, 4, 5, 6], agent-based sensus on the validity of the imbalance of regula- models [7, 8, 9, 10] and stochastic dynamics tory T cells as a therapeutic target or tool for [11, 12, 13, 14], to statistical mechanical models prognostic monitoring has not been achieved. [15,16,17,18,19,20],networkapproaches[21,22] and machine learning [23, 24, 25, 26, 27]. For re- While many mathematical models have been cent reviews, see [28, 29]. One of the reasons, is developed to understand how regulatory T cells that the rapid improvement of data storage ca- promote immunological tolerance towards self- pacity and experimental techniques, has led to an antigens (see[42,43]forareview), anunderstand- explosion of data, meaning that quantitative tools ingoftheeffectoftheregulatoryTcellsimbalance are needed to interpret and efficiently extract in- on immunosuppressionhaslagged behinddataac- 1 quisition. In this work we propose a mathemati- the germinal centre1, where B cells proliferate at cal model of the adaptive immune system that is a rate that is unparalleled in mammalian tissues. aimed at filling this gap. In Section 2 we summa- Random mutations during clonal expansion cause rizethekeymechanismsthatregulatetheadaptive the production of antibodies with a broad range immune reponse. In Section 3 we set-up a mini- of binding affinities for their antigen. B cells with malmathematical modelfor theadaptive immune unfavourablemutationswillnotgetsufficientlyac- system, comprisingconventional and regulatory T tivated by the antigen and T cells, and will die, cells and B cells and we study its response to an while those with improved affinity will be stim- antigen, in a wide range of the parameters space ulated to clone themselves. This leads to an ef- and for different mechanisms for regulatory T cell fectiveaffinity-dependentselectionprocess,within activation. Finally in Section 4 we summarize our the germinal centre, known as affinity maturation results and propose pathways for future work. [44],whichresultsinanincreaseofaffinity, bysev- eral orders of magnitude, over a time of one week in mammalians. 2 The adaptive immune system Tregcellsturnoffantibodyproductionandsup- press the immune response. The mechanistic de- The adaptive immune system is a complex net- tails of Treg cells functioning are still debated work of cells that work together to defend the [45, 46]. For example, it is not clear whether their body against pathogens such as bacteria, virus or activation relies on antigens, in the same fashion tumor cells. The most important types of cells as conventional Th cells, and whether Treg cells that constitute it, are the lymphocytes, which can can directly suppressBcells or whetherthey must be divided into B and T cells. suppressThcellsinordertosuppressBcells. Sup- B cells are equipped with B cell receptors pression of Th cells has been extensively modelled (BCR), which are able to bind to specific in the context of tolerance (see e.g. references in pathogens. WhenaBcellcomesacrossapathogen [42, 29]). In this work we model the direct sup- that can bind its receptors, it ingests it, par- pression of B cells, which has been suggested in a tially degrades it, and exports fragments of it, i.e. numberof recent studies (see e.g. [47,48, 49,50]). antigens, to the cell surface, where they are pre- sentedinassociationwithproteinsknownasMHC molecules. Each B cell has a number n = (104) 2.1 The main reactions O of identical BCR on its surface. TcellscanbedividedintoTkillers,thatdestroy As explained above, theactivation of theadaptive infected or cancer cells, and T helpers, who have immunesystem involves several reactions between theroleofsignalingBcellsandinitiateanimmune lymphocytes and other immune cells, mediated or response. Similarly to B cells, T helper cells are not by the antigen. Here we summarize the reac- equippedwithreceptors(TCR),butincontrastto tions that we shall include in our model: B cells, they do not bind to the antigen directly, but to the MHC-antigen complex on the surface Antigens (Ag) replicate at a rate r • of antigen presenting cells (APC), such as B cells and other immune cells (e.g. dendritic cells). T Ag r 2Ag (1) → helper cells have been commonly divided into T regulatory (Treg) cells and conventional T helper B-cells bind to antigens via B cell receptors (Th) cells. When a Th cell binds an APC, it gets • to form antigen presenting cells at a rate π+ activated and starts proliferating while releasing which depends on the affinity between the proteins called cytokines, which activate B cells. BCR and the antigen Once activated, B cells undergo clonal expan- sion, i.e. form many copies of identical cells shar- Ag+B π+ APC (2) ing the same antigen receptors, and secrete anti- → bodies, i.e. a free form of those receptors, that 1Germinal centres are transient structures that form can recognize and neutralize the antigen. Clonal within peripheral lymphoid organs in response to T cell- expansion usually involves migration of B cells to dependent activation. 2 Tcells get activated whenthey bindanAPC, each labelled by µ = 1,...,P. We assume that • at a rate W each B clone µ is able to bind to one type of anti- gen (for simplicity assumed to have a single epi- T +APC W T⋆+APC (3) tope) also labelled by µ. We denote by ψ , b → µ µ and p , the population densities of antigen type µ APCs are eliminated at a rate π− (e.g. by µ, B clone µ and B cells presenting antigen µ (i.e. • macrophages or via apoptosis) APC of type µ), respectively. Experimental lym- APC π− 0 (4) phocytes counts, estimate the total number of B → cells to be αN with α 1/2. We will assume that ≃ the number of B clones is sub-extensive in N, i.e. activatedTcellsstimulateBclonalexpansion • P = αNγ, with γ (0,1), so that each clone will ∈ T⋆+B λ T⋆+2B (5) contain on average a number (N1−γ) of B cells.3 O → It is then convenient to define clonal densities as or repress it b = B /cN1−γ where B is the size of B clone µ µ µ µ and c = (N0) is a constant that will be de- T⋆+B λ T⋆ (6) fined later, aOnd similarly for p = P /cN1−γ and → µ µ ψ = Ψ /cN1−γ. We can then write the following depending on the nature (Treg or Th) of the µ µ equations, describing reactions (1), (2) and (4) T cell d Bcellsarekeptinastateofactivated apopto- p = π+ψ b π−p (8) • dt µ µ µ µ− µ µ sis while undergoing clonal expansion in the d germinal centre and compete for survival sig- ψ = r ψ π+ψ b (9) dt µ µ µ− µ µ µ nal from T helper cells (germinal centre se- lection) For each T cell i, we introduce a variable ηi which B +B δ 0. (7) takes values 1 if i is helper and 1 if it is reg- → − ulator. We assume each η to be identically and i We will neglect, in our model, the presence of independently sampled from other immune cells that can serve as APC (e.g. 1+ǫ 1 ǫ dendritic cells) and we will neglect several other P(η) = δ + − δ (10) η,1 η,−1 2 2 processes, such as antigen mutation, clonal ex- pansion of T cells, activation (by T helper cells) The parameter 1 ǫ 1 quantifies the imbal- − ≤ ≤ of T killers, macrophages and other types of im- ancebetweenTh(η = 1)andTreg(η = 1)andis − mune cells, production of antibodies by B cells, directly related to the T-helper/T-regulator lym- and affinity maturation of B cells in the germinal phocyte ratio R, measured in experiments, via centre. We will simply assume that, due to this 1+ǫ process, the affinity π+ between BCR and anti- R = . (11) 1 ǫ gen is an increasing function of time. Despite its − simplicity, we believe that this modelcaptures the For ǫ = 0 one has that Th and Treg are present in basic principles of an immune responseto an anti- equal proportions, both equal to 1/2 (i.e. R = 1), genand,tothebestofourknowledge,itisthefirst forǫ > 0Thcells aremorethanTregcells (R > 1) quantitative model that reveals the ratio between and for ǫ < 0 Th cells are less than Treg cells Th and Treg cells as the single most important (R < 1). parameter that affects its behavior. We representthestate ofeach Tcell i(active or inactive) with a variable σ which takes value 1 if i i is active and 0 otherwise. Helper T cells get ac- 3 The mathematical set-up tivated via reaction (3), i.e. when their receptors We consider a population of T cells, each labelled 3The proposed scaling is consistent with experimental by i = 1,...,N, and a population of B clones2, countsoflymphocytes,thatestimatethenumberofTlym- phocytestobeapproximatelyN =O(1012),thenumberof 2A B clone is the ensemble of all B cells that have the B clones to be P =O(108), and the size of a clone at rest same receptors and thusrespond tothe same antigen. to be O(104). 3 B cell q ✗✔B clone bind to an APC. Assuming for simplicity that T T cell cells updatetheir state at regular time intervals of ✈ ✖✕ duration ∆, we have µ 1 ... µ ... αNγ σ (t+∆)= θ( ζ p (t) Tz(t)) (12) i i µ − µ ✗✔✗✔✗✔✗✔✗✔ X qq q q q q q q q wθindwTs(oihhsxraiteesm)rtrzihetab=ehelruierθost0(-ei(nxaod1ofvn))oievsrrpaoiesa(rrxgzsiantte)≤nrh,doecetwner0g(,has0tntidahc)edrnhp,aTodwammfncnθuidem(vnlfxalrcζi)orctiµiiikmao=cs∈bnaslb,ne{1oi1motbwf,alooei0kinrgt}isdhinycxigamnaslu>dmcviintecaeoal0alltuib,tsreieleinyzcss, ✡✄(cid:0)✡✄✄(cid:0)✄✡✄✄✡✄✖(cid:0)✄✄✡q(cid:0)✡✡❈(cid:0)❈✕✡❈❈❈(cid:0)✡❈❈❈✡❈✖(cid:0)☞✜☞✜☞(cid:0)☞q❉❉✜☞❉(cid:0)❉☞❉✕❉✜☞❉❉☞(cid:0)❉❉☎☎✜☞☎☎☞☎✖☎ζ✜☎☎iµ☎☎q✜q▲✜✡▲✄✄✡✕▲✜✄✄▲✡✄▲✄✜✄▲✡✄✖✄▲✄✡▲▲▲✡q☞☞✡☞☞✡❉✕☞❉❉☞❉✡❉☞❉❉✡☞❉❉✖❉☞☞q❉❉❉❉❉✕❉❉❉❉❉ µ ✈ ✈ ✈ ✈ ✈ ✈ ✈ ✈ B clone µ. We assume the variables ζ to be { i } 1 ... i ... N identically and independently distributed accord- ing to Figure 1: Schematic representation of the inter- µ c c actions between T cells and B clones. A link ζ p(ζiµ)= Nγδζiµ,1+ 1− Nγ δζiµ,0 (13) between T cell i and B clone µ means that T ceill (cid:18) (cid:19) i can bind to B cells in clone µ. The number of T so that each T cell i can signal an average num- cells signaling a clone µ is of thesame orderas the ber P ζµ = αc of B clones and each B clone h µ=1 i i number of B cells in clone µ, assumed (N1−γ). µ can be signaled, on average, by N ζµ = O P h i=1 i i In the figure, these numbers are chosen as (1) cN1−γ T cells. Since the number of B cells in O and identical and are meant for illustration only. P each clone is on average (N1−γ), the choice (13) O corresponds to the scenario where the number of n−1 of the third term accounts for the fact that T cells that can signal a B clone is of the same upon binding an antigen, a B cell only engages order as the number of B cells in that clone. Bi- one of its n receptors, so that the resulting APC ologically, this would seem the most efficient and is still a Bcell, with one sparereceptor less, which economical scenario, asitavoids, ontheonehand, leads to a decrease of the effective number of B a redundancy of B cells that would not get suffi- cells by a fraction 1/n only. In terms of popula- ciently signaled by a comparatively small number tion density, we have of compatible T cells, and on the other hand, a waste of T cells in signaling a comparatively small db 1 δ π+ µ µ µ µ number of B cells. For a schematic representation =λ b ζ η σ b ψ of the signalling between T cells and B clones see dt µ µ cN1−γ i i i i − λµ µ− λµn µ! X Figure 1. (15) Finally, B clones expand (contract) when they where the growth term is now proportional to the receive excitatory (inhibitory) signals by active T density of the net excitatory signal received by B cellsandcompeteforsurvival,sothateachBclone clone µ from T cells, that we shall denote as follows a logistic dynamics 1 dBµ = Bµ λ ζµη σ δ B 1π+Ψ mµ(σ) = cN1−γ σiηiζiµ (16) dt cN1−γ µ i i i− µ µ−n µ µ! Xi i X (14) withσ = (σ1,...,σN)representingthemicrostate (active or inactive) of all T cells. Finally, bearing wherethedenominator in the prefactor represents in mind that n = (104), and neglecting (n−1) O O the carrying capacity of the system for clone µ in (15), we get whenallcognateT-helpersareactivatedandsend- ing excitatory signal to clone µ and Ψ represents d δ µ b = λ b m(σ) µb (17) the number of antigens of type µ. The prefactor dt µ µ µ − λµ µ! 4 Denoting by (z x) = x dzp(z) the cumula- W (σ ) at which the cell is flipped. Multiplying P ≤ −∞ t j µ tive distribution function of the noise distribution (19) times η ζ and summing over j, we obtain R j j p(z), it follows from equation (12) that the likeli- the following equation of motion for the average hood to observe configuration σ at time t+∆, is, signal strength m (t) = m (σ(t)) on clone µ: i µ µ h i for any symmetric distribution p(z) = p( z), − dm Nγ β µ = ζµη[1+tanh ζνp ] m dt 2c h 2 ν iη,ζ − µ µ ν pt+∆(σi)= P z ≤ (2σi −1)β µ ζi pµ(t)! X (20) X In the above ... denotes the average h iη,ζ where β = 1/T is the inverse noise level. A natu- ...P(η,ζ)overthedistributionofregulatory η,ζ ral choice for p(z) would be a Gaussian distribu- patterns in the system (46). P tion, that leads to (z x) = 1(1+erf(x/√2)). Assuming that P(η,ζ) = P(η)P(ζ) i.e. the P ≤ 2 An alternative choice is the so-called Glauber dis- ability of a helper cell i to bind to clone µ does tribution (commonly made in statistical physics) not depend on whether i is a helper or regulator, leading to we have 1 x dm ǫNγ β (z x) = 1+tanh , µ = ζµ[1+tanh ζνpν] ζ mµ. (21) P ≤ 2 2 dt 2c h 2 i − (cid:18) (cid:19) ν X qualitatively very similar to the cumulative distri- Finally, averaging over ζµ,usingtheindependence bution function for the Gaussian distribution and of the ζµ’s and ζµ = c/Nγ, we get h i easier to work with analytically. For this choice, the probability that helper T cell i changes in a dm ǫ β µ= 1+ tanh p + ζνp m , single time step its state σi at time t (to 1 σi at dt 2 h 2 µ ν i{ζν}− µ time t+∆) is − (cid:16) νX6=µ (cid:17) (22) where one has still to carry out the average over 1 β Wt(σi) = 2"1+(1−2σi)tanh 2 µ ζiµpµ(t)# tnhaeiv{eζlyν}exopthecetr tthhaatneζqµu.atAiotnfi(r2st2)gldaonecse,noonteclmosaey, X as m depends on the whole set of p , and each where we used 1 2σi = 1 and tanh( x) = µ { ν} − ± ± p depends, via (8), on b and ψ . The latter is tanhx. Assuming that the update of T cells is µ µ µ ± determined from b , however the former depends, sequential, i.e. at each time step one helper cell µ via(17),onm (σ),whichisaprioridifferentfrom i, drawn at random, is updated with likelihood µ the m sitting in (22), defined as the average of W (σ ), one obtains, for ∆ = 1/N and N large, µ t i m (σ). However, one can show (see Appendix A) the following master equation for the probability µ that for γ < 1, fluctuations of m (σ) about its density p (σ) to observe microstate σ at time t, µ t thermodynamic average m vanish for large N. µ d This allows us to replace m (σ) in (17) with m p (σ)= [p (F σ)W (1 σ ) p (σ)W (σ )] µ µ t t i t i t t i dt − − giving i X (18) db δ where Fi is a “cell-flip” operator that changes the µ = λµbµ mµ µbµ (23) configurationofTcellifromσ to1 σ andhasno dt − λµ ! i i − effect on any other cell j = i. From (18) one can This finally enable us to express the dynamical 6 derive equations of motion for the expectations evolution of this immune system model in terms = σ pt(σ). Multiplying (18) by σj and of a closed set of first order differential equations, h · i · summingover σ, weobtain thefollowing equation P namely (8), (9), (22) and (23). Equation (8) im- d plies that the equilibrium density of antigen pre- σ = (1 2σ )W (σ ) (19) j j t j senting cells is given, for all µ, by dth i h − i intuitively meaning that the rate of change of the p = a b ψ (24) µ µ µ µ average activity of Tcell j is given by its variation where a = π+/π− quantifies the affinity between 1 2σ upon a single cell flip F , times the rate µ µ µ j j − 5 0.5 b ψ µ 0.4 µ 800 0.3 600 0.2 400 0.1 200 0.0 0 10 20 30 40 50 10 20 30 40 50 t t (cid:0)0(cid:4) 1000 (cid:0)0(cid:3) 800 (cid:0)0(cid:2) 600 (cid:0)0(cid:1) 400 (cid:0)01 200 (cid:0)0(cid:0)0 20 30 40 50 10 20 30 40 50 1(cid:0) Figure 2: Time evolution of B clonal density b (left) and antigen concentration ψ (right) for r = 2, µ µ µ ǫ = 0.5, λ = δ = π− = 1 and π+ = 10 [tanh(10)+tanh(t 5)]. Top panels: β = 10. Bottom µ µ µ µ × − panels: β = 0.1; Initial conditions were chosen as ψ (0) = 0.3, b (0) = κ m (0), m (0) = ǫ/2 and µ µ µ µ µ p (0) = 0. For higher noise levels T = β−1, the system mounts a weaker immune response, resulting µ in lower B clonal densities and higher antigen concentrations, however the antigen is still removed within the same timescale. Bclone µ and antigen µ.4 If thesystem is initially as for ψ = 0, p = 0 is a fixed point of the dy- ν ν in equilibrium in the absence of antigen, one has namics (8), and from (22) we get p = 0 µ, and from (22) we get m = ǫ/2 µ. µ µ dm ǫ β ∀ ∀ ν = 1+ tanh ζµp m Hence, equation(23)yieldsbµ = κµǫ/2, withκµ = dt 2 h 2 µiξµ − ν (cid:20) (cid:21) λµ/δµ, for ǫ > 0, and bµ = 0 for ǫ < 0. = ǫ m + (N−γ), (26) ν 2 − O 3.1 Response to a single antigen of which m = ǫ/2 is a fixed point, for large N. ν In contrast, equation (22) gives for clone µ, that In this section, we study the activation of the is compatible with the antigen introduced in the immune system when a single antigen type µ is host, present, i.e. ψ = 0 and ψ = 0 ν = µ. µ ν 6 ∀ 6 Assumingthat, when antigen µ is introduced in dmµ ǫ β = 1+tanh p m (27) µ µ the host at time t = 0, the latter is resting at dt 2 2 − (cid:20) (cid:21) its equilibrium state in the absence of antigen, we whose stationary solution is, inserting (24), have for all t > 0 and all ν = µ ǫ β 6 m = 1+tanh a b ψ (28) ǫ κ ǫ µ µ µ µ ν 2 2 (ψν,pν,mν,bν) = 0,0, ,max 0, , (25) (cid:20) (cid:21) 2 2 (cid:18) (cid:26) (cid:27)(cid:19) where bµ and ψµ are the steady state solutions of 4Theaffinitycanbeexperimentallymeasuredasthera- (23) and (9), respectively. Equation (28) shows tio between theequilibrium concentration of APC and the that at stationarity m > 0 for any ǫ > 0, so the µ product of the equilibrium concentrations of antigen and stable steady state of (23) is b = κ m . specific B cells. µ µ µ 6 0.5 b ψ 80 µ 0.4 µ 60 0.3 40 0.2 0.1 20 0.0 0 10 20 30 40 10 20 30 40 t t 0.5 400 0.4 0.3 300 0.2 200 0.1 100 0.0 0 10 20 30 40 10 20 30 40 Figure 3: Time evolution of B clonal density b (left) and antigen concentration ψ (right) for r = 1, µ µ µ ǫ = 0.5, λ = δ = π− = 1 and β = 1. Top panels: π+ = 10 [tanh(10) +tanh(t 7)]. Bottom µ µ µ µ × − panels: π+ = 10 [tanh(2.01)+tanh[(t 20)/10]]. Initial conditions were chosen as in Figure 2. The µ × − time-dependence of the affinity maturation π+ affects both the intensity of viral concentration and µ the timescale on which viral removal is accomplished. This yields, for the stationary value of m , the for π+ < r /ǫκ , the viral concentration will grow µ µ µ µ solution of the self-consistency equation indefinitely, whereas for r /ǫκ < π+ < 2r /ǫκ , µ µ µ µ µ the system will settle, after a transient, onto a ǫ βκ mµ = 1+tanh µaµψµmµ (29) limit cycle, where mµ(t) oscillates between two 2 2 (cid:20) (cid:21) values m−, m+, with ǫ/2 m− < m+ ǫ, such ≤ ≤ where ψµ is the steady state solution of that rµ − κµπµ+m− > 0 and rµ − κµπµ+m+ < 0. This leads to a periodic motion of viral concen- d ψ = ψ r κ π+m . (30) tration, meaning that the antigen will permane dt µ µ µ− µ µ µ indefinitely in the host, but its growth is limited (cid:16) (cid:17) by the action of the immune system. One can see that (ψ ,m ) = (0,ǫ/2) is the only µ µ stable fixed point of this dynamical system for As B cells undergo clonal expansion, they in- π+ > 2r /ǫκ , meaning that for ǫ > 0, anti- crease their affinity with the antigen, via affin- µ µ µ gen removal is accomplished when the affinity of ity maturation, so πµ+ is an increasing function the antigen-specific BCR is large enough com- of time. For the model in study here, this implies pared to the antigen replication rate. Conversely, that, as long as ǫ > 0, the system will eventu- for πµ+ < 2rµ/ǫκµ, the fixed point (0,ǫ/2) be- ally reach the affinity value πµ+ = 2rµ/ǫκµ where comes unstable and different types of dynamics the fixed point (ψµ,mµ) = (0,ǫ/2) becomes stable may arise, depending on the range of the kinetic and thesystem will manage toremove theantigen parameters. Since 0 tanh(βa b ψ /2) 1, one completely. µ µ µ ≤ ≤ has, from (28), ǫ/2 m ǫ, hence the largest In contrast, for ǫ 0, we get from (28) m 0, µ µ ≤ ≤ ≤ ≤ valueattainablebym isǫ. Then(30)impliesthat so the stable fixed point of (23) is b = 0, which µ µ 7 5 80 b ψ µ 4 µ 60 3 40 2 20 1 0 10 20 30 40 10 20 30 40 t t 5 80 4 3 60 2 40 1 20 0 10 20 30 40 10 20 30 40 Figure 4: Time evolution of B clonal density b (left) and antigen concentration ψ (right) for β = 1, µ µ π+ = 10 [tanh(10)+tanh(t 5)], π− = 1, r = 2 and ǫ = 0.5. Top panels: λ = 10 and δ = 1. µ × − µ µ µ Bottom panels: λ = 1, δ = 0.1. Initial conditions were chosen as in Figure 2 µ µ inserted into (28) and (9), leads to m = ǫ /2 that the ratio π+( )/π+(0) = (π + 1)/(π µ −| | µ ∞ µ 0 0 − and to an antigen indefinitely replicating in the tanhγt⋆) between the final and the initial value host through dψ /dt = r ψ , no matter how fast of π+ is (104 105). However, the actual time- µ µ µ µ O − or large π+ is growing. dependence of π+ has no qualitative impact on µ µ We show the results of numerical solution of the immune response, as long as it is an increas- (23), (9), (27) and (8) in Figures 2, 3, 4 and 5, ing function, whose limit is π+( ) > 2r /ǫκ . µ ∞ µ µ where b and ψ are plotted as a function of time Plots show an immune response that follows per- µ µ for different choices of the inverse noise level β fectly the behaviour of the antigen, i.e. it ex- (Fig. 2), time-dependence of π+ (Fig. 3), differ- pands B clones while the antigen concentration µ ent values of the kinetic parameters λ , δ (Fig. is increasing and contracts them when the anti- µ µ 4),andfordifferentvaluesofπ− andantigenrepli- gen concentration is decreasing, so that both the µ cation rate r (Fig. 5). On the assumption that B clone concentration and the antigen concentra- µ the affinity is an increasing function of time that tion are unimodal functions of time. The stochas- saturates at large time, π+(t) was chosen as the tic noise T has the effect of mildly reducing the µ following sigmoid function: height of the peak in B clones concentration, thus increasing the peak in viral concentration, how- π+(t) = π [π +tanh(γ(t t⋆))] (31) µ M 0 − ever the system will be able to remove the antigen where γ,t⋆,π and π are parameters that con- even at high noise levels (see Fig. 2). The time- M 0 trol, respectively, how fast, early and large the dependence of πµ+, affects both the location and affinity grows and its initial value. Since the theheightofthepeaks,consistently withtheintu- affinity is experimentally found to increase by ition that the faster the affinity grows, the earlier four or five orders of magnitude over a pe- and the smaller the peak in viral concentration. riod of a week, we chose t⋆ 7 and π such Fig. 3 shows that the affinity that starts increas- 0 ≤ 8 0.4 1.5 b ψ µ µ 0.3 1.0 0.2 0.5 0.1 0.0 0 10 20 30 40 50 10 20 30 40 50 t t 0.5 2.0×109 0.4 1.5×109 0.3 1.0×109 0.2 5.0×108 0.1 0.0 0 10 20 30 40 50 10 20 30 40 50 Figure 5: Time evolution of B clonal density b (left) and antigen concentration ψ (right) for β = 1, µ µ π+ = 10 [tanh(10)+tanh(t 5)], λ = δ = 1 and ǫ = 0.5. Top panels: r = 0.5 and π− = 0.1. µ × − µ µ µ Bottom panels: r = 5, π− = 1. Initial conditions were chosen as in Figure 2. µ ing at day t⋆ = 7 with rate γ = 1 outperforms More in general, this model captures the effec- the affinity that starts increasing at day t⋆ = 20 tiveness of the adaptive immune system in deal- with a slower rate γ = 1/10. The kinetic parame- ing with very diverse scenarios and predicts that ters λ and δ have a mild effect on the shape of the most important single parameter in determin- µ µ the B clonal concentration (see Fig. 4), while Fig. ing whether the system is in a healthy or in an 5 (top panel) shows that π− affects more notice- immuno-suppressed phase is ǫ, directly related to µ ably the decay of B cells after the infection peak, the Th/Treg ratio R via (11) and in line with re- potentially sustaining long-term memory. Finally, centexperiments[36,37,41]. Finally,wenotethat Fig. 5 (bottom panel) shows that although faster- althoughunstable,b =0isalwaysafixedpointof µ replicating antigens will attain much higher con- (23), even for ǫ > 0. This suggests that a fast mu- centrations,theyarestillremoved,forpositiveval- tating pathogen (like e.g. cancerous cells or HIV) ues of ǫ, within similar timescales, by triggering a that manages to release new (i.e. mutated) types stronger immune response. of antigens in thehost, beforebeingremoved from In contrast, Fig. 6 shows that as soon as ǫ the system, may succeed in escaping the immune is lowered below zero, the system becomes un- system, by eventually mutating to a form µ⋆ for responsive and is not able to fight a single anti- which there is no answering B cell, i.e. bµ⋆ = 0. gen, even if replicating slowly. Ourresults suggest Given that B cells responding to self-antigens are that for Th/Treg ratios R > 1 the host manages suppressed via negative selection processes occur- to remove completely the antigen and the time it ring in the bone marrow, one way for a non-self takes depends on how large and fast the affinity antigen to escape theimmune system is to mutate π+ between BCR and antigen grows. Conversely, to a form close to self-antigens, as it occurs in µ for R < 1 the immune system is impaired and HIV. Another way is to disruptthe mechanism by infections will permane indefinitely in the host. which non-self antigens are flagged up as such, by 9 0.10 b ψ µ 0.08 µ 3 0.06 2 0.04 1 0.02 20 40 60 80 100 20 40 60 80 100 t t 0.10 1.2×1018 0.08 1.0×1018 0.06 8.0×1017 6.0×1017 0.04 4.0×1017 0.02 2.0×1017 20 40 60 80 100 20 40 60 80 100 Figure 6: Time evolution of B clonal density b (left) and antigen concentration ψ for β = 1, µ µ π+ = 10 [tanh(10)+tanh(t 5)], π− = 1, λ = 5, δ = 1, r = 0.5. In the top panels ǫ = 0.01, µ × − µ µ µ µ while in bottom panels ǫ = 0.01. Initial conditions were chosen as ψ(0) = 0.3 and b (0) = κ ǫ /2. µ µ − | | The plot shows that as soon as ǫ < 0 the system cannot remove the antigen. e.g. altering the functionality of MHC molecules, by a random variable, while for η = 1 we have i a pattern that has been observed both in cancer the same equation as before. Repeating the steps and HIV progression [52, 53]. This secound route that led to (22), we get µ would make ζi = 0 ∀ i, so that mµ = 0 and the dmµ = Nγ ζµη[1+tanh β1+η ζνp ] ] m density of responding B cells becomes bµ = 0. dt 2c h 2 2 ν iη,ζ − µ ν X (33) 3.2 The role played by the mechanism that simplifies, in the single antigen case, to that activates Tregs dm 1 1+ǫ β µ = ǫ+ tanh p m (34) µ µ The activation of Treg cells is still quite debated. dt 2 2 2 − (cid:20) (cid:21) Above we assumed that Treg cells get activated For ǫ > 0, the stationary solution of (34) is in the same way as Th cells, i.e. by binding to m > 0, hence the stationary solution of (23) is µ an APC. However, there is no general consensus b = κ m where m solves the self-consistency µ µ µ µ about the need of antigen mediation [54]. equation This prompts us to consider an alternative ver- 1 1+ǫ β sion of the model used above, where the mech- mµ = ǫ+ tanh aµψµκµmµ (35) 2 2 2 anism that activates Th cells remains the same, (cid:20) (cid:21) Thesolution is easily foundin terms of theshifted while activation of Treg cells is fully stochastic variablesM = m ǫ/2byintersectingthecurves µ µ − 1+η 1+ǫ i µ σ (t+∆) = θ ζ p (t) Tz(t) (32) M = tanhx (36) i 2 i µ − ! µ 4 µ X 2x ǫ M = (37) i.e. for ηi = −1 the activation of i is determined µ βκµaµψµ − 2 10