Physiological factors leading to a successful vaccination: A computational approach

Journal of Theoretical Biology 454 (2018) 215–230

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/jtbi

Physiological factors leading to a successful vaccination: A computational approach Dominic L. Maderazo a, Jennifer A. Flegg a, Melanie R. Neeland b, Michael J. de Veer c, Mark B. Flegg d,∗ a

School of Mathematics and Statistics, The University of Melbourne, Australia Murdoch Children’s Research Institute, Australia Department of Physiology, Faculty of Medicine, Nursing and Health Sciences, Monash University, Australia d School of Mathematical Sciences, Monash University, Australia b c

a r t i c l e

i n f o

Article history: Received 16 March 2018 Revised 1 June 2018 Accepted 6 June 2018

Keywords: Immune system Multiscale mathematical modelling Computational simulation Lymph node

a b s t r a c t The immune system mounts a response to an infection by activating T cells. T cell activation occurs when dendritic cells, which have already interacted with the pathogen, scan a T cell that is cognate for (responsive to) the pathogen. This often occurs inside lymph nodes. The time it takes for this scanning event to occur, indeed the probability that it will occur at all, depends on many factors, including the rate that T cells and dendritic cells enter and leave the lymph node as well as the geometry of the lymph node and of course other cellular and molecular parameters. In this paper, we develop a hybrid stochasticdeterministic mathematical model at the tissue scale of the lymph node and simulate dendritic cells and cognate T cells to investigate the most important physiological factors leading to a successful and timely immune response after a vaccination. We use an agent-based model to describe the small population of cognate naive T cells and a partial differential equation description for the concentration of mature dendritic cells. We estimate the model parameters based on the known literature and measurements previously taken in our lab. We perform a parameter sensitivity analysis to quantify the sensitivity of the model results to the parameters. The results show that increasing T cell inflow through high endothelial venules, restricting cellular egress via the efferent lymph and increasing the total dendritic cell count by improving vaccinations are the among the most important physiological factors leading to an improved immune response. We also find that increasing the physical size of lymph nodes improves the overall likelihood that an immune response will take place but has a fairly weak effect on the response rate. The nature of dendritic cell trafficking through the LN (either passive or active transport) seems to have little effect on the overall immune response except if a change in overall egress time is observed. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction The immune response is the body’s biological mechanism to defend against infections and comprises of both innate and adaptive immunity. While innate immunity, the first line of defence, is often sufficient to control an infection, developing essential long-term immunity requires an adaptive response. The adaptive immune response is multifaceted, however, central to the process is the activation of immune cells known as T cells. Inactivated (or naive) T cells patrol through blood and lymphatic systems and are activated by the detection of antigen, a marker for a specific infection. T cells express receptors (TCRs)

∗

Corresponding author. E-mail address: mark.fl[email protected] (M.B. Flegg).

https://doi.org/10.1016/j.jtbi.2018.06.008 0022-5193/© 2018 Elsevier Ltd. All rights reserved.

which are only specific to an associated (or cognate) antigen. T cell activation occurs when TCRs engage with a processed form of their cognate antigen known as pMHC (antigenic peptide bound to major histocompatibility complex) (Kindt et al., 2007). Efficient activation of T cells requires macrophages and dendritic cells (DCs) which express large numbers of MHC proteins. Activation causes rapid T cell proliferation and the daughter cells subsequently mount an antigen-specific response (Kindt et al., 2007). Each T cell expresses a specific TCR and each organism needs many T cells and a wide variety of TCRs to help identify a wide range of pathogenic threats. The organism therefore has a TCR repertoire. Talking about T cells, TCRs, cognate T cells/TCRs can get confusing without a well defined language to describe them. We will use the terms defined by Laydon, Bangham and Asquith (Asquith et al., 2015). The total number of distinct TCRs which

216

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

are created by an organism we will refer to as the TCR ‘diversity’ and we will use the term ‘species’ to describe the collection of T cells which share a common TCR (also known as clonotypes). We shall use the symbol NT CR to describe the TCR diversity and the symbol Nt to describe the copy number of clonotypes of a given species. The copy number of clonotypes of the various species of T cell in an organism is heavily skewed but also specific to the individual. We treat Nt therefore as the effective ‘average’ copy number of clonotypes for each species. We can relate Nt and NT CR therefore with the total number of T cells in the whole organism Ntotal = NT CR Nt . As an example, for the αβ T cell repertoire in mice the diversity is estimated at NT CR ∼ 2 × 106 (Casrouge et al., 20 0 0) whilst humans have genes which encode for a larger number of TCRs NT CR ∼ 2 × 107 (Naylor et al., 2005). Estimates for the total number of naive T cells in humans vary depending on the body mass and age of the individual, but are on the order of Ntotal ∼ 1011 (Bains et al., 2009) and therefore a given T cell species has approximately Nt ∼ 10, 0 0 0 clonotypes. Finding an antigen-specific naive T cell (a T cell of a specific species) amongst the total pool of T cells is like finding a needle in a haystack. T cell activation is supported by dendritic cells (DCs). After processing antigen at the site of an infection, DCs undergo a maturation phase during which they grow long dendrites which can extend up to twice their body length (Miller et al., 2004b). Fragments of the processed antigen are displayed on MHC molecules found on the outer membrane of the DCs and the purpose of the dendrites is to interact with as many neighbouring T cells as possible, optimising the likelihood that a T cell specific to the antigen will be primed (von Andrian and Mempel, 2003; Randolph et al., 2005). For this reason, DCs are often called antigen presenting cells (APCs). Lymph nodes (LNs) are secondary lymphoid tissues. Their main function is to facilitate the otherwise unlikely event that an antigen-specific T cell (a T cell belonging to a species cognate for the exposed antigen) will encounter its cognate antigen on a mature DC and be activated. Until the mid 20 0 0s, theoretical studies of LNs considered them to be spatially homogeneous in which naive T cells could encounter mature DCs in a closed, concentrated environment at a rate determined only by copy number (and not spatial distribution). However, LNs are complex tissues with a spatially inhomogeneous internal structure that supports the interaction between DCs and T cells, although how this is achieved is poorly understood. A schematic cross-section of the anatomy of a LN is presented in Fig. 1. Our schematic is based on detailed diagrams presented by Willard-Mack (2006). At the site of an infection (or vaccination), immature DCs ingest antigen, mature and are actively transported through lymphatic vessels to LNs. The lymph (the fluid in the lymphatic vessels), which carries the mature DCs, drains into a nearby LN and enters the capsule (‘top’) of the LN via one of the afferent lymphatic (AL) vessels (Fig. 1). Lymph that enters the LN flows into the subcapsular sinus and preferentially around the lobules, through transverse sinuses, towards the efferent lymphatic (EL) vessel where it leaves the LN. The DCs which are carried into the LN via the AL deposit on the interface between the superficial cortex and the subcapsular sinus. Scanning electron micrographs of rat LNs have revealed that the floor of the subcapsular sinus possesses pores leading to the paracortex (Ohtani and Ohtani, 2008). DCs are transported rather directly through the superficial cortex, which consists largely of follicles filled with B cells (which are not the focus of this manuscript), into the paracortex and the deep cortical unit (DCU) which contains densely packed naive T cells which are constantly being trafficked through the LNs of the organism. Twophoton experiments have provided evidence that the DCs, at a cellular level, gradually move through the cortex, in a random fashion (Grigorova et al., 2010; Miller et al., 2002). This microscopic move-

Fig. 1. Schematic cross-section of a lymph node (LN). Lymph drains DCs into the LN from nearby tissues via an afferent lymphatic (AL) vessel and then flows through the subcapsular and transverse sinuses to the efferent lymphatic (EL) vessel, whilst depositing DCs into the lobules. Lobules consist of two main regions; the superficial cortex (between green and yellow lines) typically consists of several follicles which contain B cells (maroon) and the paracortex (between yellow and orange lines), often referred to as the T cell zone because it contains mostly T cells. Within the paracortex, the deep cortical unit (DCU) contains a concentrated population of T cells (dark blue); it is in the DCU that the majority of TC-DC interactions occur. Cells from the lobules enter the medullary sinuses and are drained into the medulla (within the orange line). T cells in the medulla drain into the EL for recirculation whilst DCs are removed by macrophages. Blood supply brings T cells into the LN, where T cells are deposited into the lobules through high endothelial venules (HEVs), shown schematically in thin red lines. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ment is consistent with Brownian motion and diffusion constants can be calculated from experimental data (Miller et al., 2004a). However, spread throughout the DCU are sinuses carrying lymph which flows from the AL to the medulla and eventually into the EL. DCs can find these sinuses and actively move with a bias towards the EL. Therefore, it is expected that there may be some net active transport of cells from the AL to the medulla. The network of sinuses may also contribute to an increase in tissue scale dispersion of cells as well as active transport from AL to EL, but it is unclear to what extent this occurs (Murphy et al., 2007; SainteMarie et al., 1982). There is no evidence that mature DCs exit the LN; they are not present in efferent lymph or circulating blood and typically reside only within the LN after entering. Mature DCs in the murine LN are known to reside for an average of about 60 h (Haig et al., 1999; Neeland, 2015). It is likely DCs are triggered to apoptose (die) within the medulla and are then scavenged by resident macrophages (Tomura et al., 2014). DC fate within the LN is rather poorly understood. However, assuming DCs are removed on exiting the DCU and entering the medulla, the time spent in the LN provides us with a more reliable source of approximate transport parameters on the tissue scale than the microscopic scale transport parameters determined from 2-photon experiments. Naive T cells originate from haematopoietic stem cells in the bone marrow and develop in the thymus before being distributed around the body. T cells typically enter a LN through the blood supply. The blood vessels branch into small vessels known as high endothelial venules (HEVs) that are distributed throughout the lobules in the LN. T cells in the HEVs adhere to the endothelial walls leading to transendothelial migration into the LN and particularly into the DCU (Sage and Carman, 2009). Interestingly, it has recently been shown that mature DCs can activate HEVs so that they become more permeable to T cells entering the DCU (Moussion and Girard, 2011). In the DCU tissue, T cells have also been shown to

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

move randomly using two photon experiments but substantially more rapidly than their DC counterparts and show no evidence of chemotactic attraction to DCs (Miller et al., 2004a). Recently, it has been shown that T cells in the DCU may enter very small lymphatic vessels known as cortical sinuses which permeate through the DCU and are often found close to HEVs (Cyster and Schwab, 2012; Förster et al., 2012). These vessels drain into the medullary sinuses where T cells are carried into the EL with the bulk flow of the lymph to recirculate through the body. The full nature of how cortical sinuses transport T cells in the DCU is not well understood. Whilst it may be assumed that T cells flow out of the LN once they enter a cortical sinus, it has been suggested that T cells may also roll along the outer surface of these small vessels before detaching later and even enter the cortical sinus for a short time before re-entering the DCU, ultimately contributing to some component of active transport of T cells in the DCU at the tissue scale in the direction of lymph flow (Cyster and Schwab, 2012; Grigorova et al., 2010; 2009). As with DCs, the transport of T cells at the tissue scale is best characterised by known macroscopic trafficking rates and dwell times (rates of inflow and outflow). There is no evidence to suggest that cognate T cells have an affinity for a particular lymph node. We will therefore consider the number of T cells inside a given DCU is some unbiased random sample of the T cells in the whole system. We will define ntotal =  Ntotal as the total (large) number of T cells in the DCU where  is the fraction of T cells present in the given DCU compared to the whole lymphatic system. We also define nt =  Nt as the expected number of T cells of a single species in the DCU at any moment in time, however as this number is usually very small we will intrinsically observe substantial fluctuations of these copy numbers. In the DCU, slowly moving DCs and relatively motile antigenspecific naive T cells can interact when a TCR from a T cell engages with the antigen presented by a DC. Despite being surrounded by naive T cells, there may only be a few T cells specific for the antigen being presented by the DC in the entire LN and the necessary encounter is subsequently unlikely. We shall use the symbol r to denote the (small) frequency/proportion of T cells which are cognate for a given antigen. Thus, NAg = rNT CR (the ‘cognate TCR diversity’) is the approximate number of distinct T cell species which are cognate for a given antigen (and for which there needs to be at least one encounter each with the cognate antigen for optimal immune response). Thus, after an infection or vaccination (which usually involves giving a patient a reduced strength pathogen to elicit an immune response) the LN typically experiences a large influx of mature DCs depending on the size/dose and type of infection/vaccination. How many mature DCs enter the LN varies greatly between infections and vaccination regimes. As DCs move slowly through the densely packed T cells they are said to ‘scan the repertoire’ of TCRs as they interact with individual T cells, that is, they scan as much of the full TCR diversity as they can to ensure activation of all cognate T cells. Naive T cells not specific to the antigen presented by a DC take approximately 3 minutes to be scanned and a single DC can scan between 500 and 5000 T cells per hour (Bousso, 2008). Combined with the fact that the DCU is highly concentrated with T cells, the purpose of the DCU in this context of the immune system is clear; the DCU substantially increases the chance that a cognate T cell will encounter a mature DC (Bousso and Robey, 2003; Miller et al., 2004a). Due to the efficiency of DC scanning capabilities, the rate of interaction between T cells and DCs is thought to be predominantly transport-limited. In order to ensure safety yet maintain immunogenicity of vaccine formulations, immunologists are often faced with the question, what is the minimal amount of antigen required in a vaccine formulation to ensure all cognate T cell species are activated. The answer to such a question depends on many factors, including transport rates of both DCs and T cells, the size of the DCU, as well as the extent

217

to which HEVs supply T cells and lymphatic sinuses remove cells from the DCU. In this paper, we will explore the effect of these factors on the immune response. A number of mathematical models have quantified DC-T cell interactions, using three main modelling approaches. The first and simplest modelling approach treats the LN as a spatially homogeneous compartment and describes the rates of change of DC and T cell concentrations using ordinary differential equations (Bajaria et al., 2002; Lee et al., 2012; Marino et al., 2004). The advantages of these models is they are simple to construct and can typically be easily solved numerically. However they do not capture the spatialheterogeneity within a LN or the stochastic nature of cognate T cell scanning events. The second modelling approach describes cells (or agents) individually and are known as agent-based models (ABM). Cells move and interact according to probabilistic rules based on the current state and position of the cell. ABMs have had many successes modelling host-pathogen systems (for a review, see Bauer et al., 2009). The main advantage of ABMs is that spatial and environmental detail within the LN can be included. Riggs et al. used an ABM to investigate the effect of HEV-to-cortical sinus average separation on scanning times and turnover times for T cells (Riggs et al., 2008) while Linderman et al. used an ABM to investigate DC scanning of T cells during acute infection on a cellular scale and inflow/outflow of cells was calibrated from human studies (Linderman et al., 2010). These sorts of models are also extensively utilized when investigating the interactions between DCs and T cells when the agents search for each other in space (Bogle and Dunbar, 2010a; 2010b; Celli et al., 2012; Gong et al., 2013). Furthermore, because of the intrinsic probabilistic nature of these models, stochastic effects are naturally incorporated into the ABM framework. A more detailed modelling approach are cellular Potts models (CPMs), which typically involve discretising the domain and individual cells occupy a number of lattice points. This means that cells have defined sizes/morphologies and rules governing motility can be made on a physical (rather than phenomenological) basis. Whilst CPMs tend to give very detailed results, the computational demand of these models are prohibitively large to model the tissue scale and instead have been used to model the interaction of subpopulations of T cells and DCs at a microscopic level (Beltman et al., 20 07a; 20 07b). CPMs have been an important tool to shed light on the dynamics of cellular motion of immune cells (MeyerHermann and Maini, 2005). A review of previous mathematical models of LN dynamics during infection is by Mirsky et al. (2011). Later ABM and CPM models focused on individual cells and were performed on scales comparable to HEV to cortical sinus separations, in contrast to earlier ODE models that treated LNs as spatially homogeneous.Furthermore, while cellular scale properties in the LN are critical to the overall function of the LN, to date there have been few attempts at a spatial model of LN function on a tissue scale. An interesting mystery that seems to not be addressed anywhere in the biological or mathematical literature is that despite the variance in the size of lymphatic organs between animals, adaptive immunity (which operates with the same fundamental cellular mechanisms) tends to occur on the same time frame (Miao et al., 2010; Neeland et al., 2014). How can a mouse and a sheep elicit an immune response on the same temporal scale despite clearly different spatial scales if T cell activation is transportlimited? Volume changes in secondary organs like the LNs can also occur within the same individual if it becomes inflamed during an immune response. Inflamed LNs commonly undergo a temporary process known as shut down which is driven by prostaglandin production within the LN and it significantly reduces the egress of cells via the EL (Hopkins et al., 1981; Neeland et al., 2017). Given

218

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

the hypothesis that phagocytosis of DCs is transport-governed, geometry changes during events such as shut down may have significant implications for DC-T cell interactions. In this paper, we construct and analyse a coarse-grained hybrid ABM model of a DCU in a LN after vaccination to see how the scanning times change due to the effect of tissue-scale morphological changes as well as cellular transport and copy number parameters in the DCU. In order to construct the model, microscopic features have been averaged and replaced by macroscopic properties. Furthermore, the model does not consider all ntotal naive T cells but rather only those that are cognate for the antigen presented by the DCs. Since this is usually a small number, we model each of these T cells individually and stochastically within the DCU. Due to the large number of DCs that flood a LN after inflammation or vaccination, we model the DCs using a partial differential equation (PDE) solved on an irregular lattice over the DCU. The manuscript is structured as follows. In the next section, the methods are detailed. This is followed by simulation results and parameter sensitivity tests, and finally, we discuss the physiological factors which play the most significant role in the immune response. It should be clear at this point that capitalised N (such as Ntotal and Nt ) are used for expected T cell copy numbers in the whole lymphatic system of the organism, calligraphic N (such as NT CR and NAg ) are used to describe a number of T cell species and lower case n (such as ntotal and nt ) as equivalent to their capitalised counterpart but specifically for the DCU rather than the whole system. 2. Materials and methods

Fig. 2. The geometry of the model DCU is a circle. The boundary ∂ a represents the interface between the DCU and the superficial cortex in the LN. Boundary ∂ e represents the interface between the DCU and medulla, where DCs are phagocytosed and T cells recirculated.

2.1. Parameters •

As definitive parameters for our model are difficult to obtain, we will focus specifically on modelling sheep lymph nodes. We have obtained cellular trafficking data for sheep lymph nodes in our lab and can be relatively certain about the magnitude of cellular copy numbers. We model only a single DCU. We do this for two reasons. First, as previously mentioned DCs are escorted rather efficiently to the paracortex and this is where the majority of TC interactions occur. It is not immediately clear that DCs spread over many DCUs after a vaccination event. Second, we are left to assume that DCUs in all of the draining LNs behave in a parallel and independent way to each other and thus we study the singular DCU structure to gain insight into the efficacy of this unit in providing an immune response. Later, we shall demonstrate that if the DCU environments are comparable, the distribution of DCs between DCUs should not affect the immune response (according to our model). In Table 1, we present parameters which are either determined from the literature, our previously published experiments or inferred later in this manuscript from our model for an ovine DCU. 2.2. Model geometry and overview Our model considers a single 2D DCU with a circular domain

 ⊂ R2 of radius R0 centred at the origin. We break the boundary ∂  into three parts, ∂  = ∂ a ∪ ∂ e ∪ ∂ ρ . In Table 1, we estimate the radius of an ovine DCU based on a full LN volume of ∼ 1.9 × 1012 μm3 from Neeland (2015) and the following assumptions/approximations: •

The large proportion of the volume of a LN is in the paracortex (made up of the various DCUs). We approximate the paracortex to be roughly the same volume as the whole LN. Here we are interested in an appropriate order of magnitude since there is substantial variation between LNs.

•

In the LN there are a ‘few’ DCUs. We assume that there are three DCUs. Therefore the volume of a DCU is ∼ 6.3 × 1011 μm3 . The DCU is roughly spherical in shape. This gives us a radius (or characteristic length scale) of R0 ∼ 5.3 × 103 μm.

Our model domain is shown overlaid on a biological diagram in Fig. 2 to provide context. Modelling in 3D is a popular choice at the cellular scale. Many 3D models focus on a detailed simulation of reality and are designed to replicate the system to high precision (Bogle and Dunbar, 2010a; 2010b; Celli et al., 2012; Gong et al., 2013). While a 3D DCU would be more physically realistic, we feel that it is unnecessary and presumes that biological parameters are accuracy known, which is most definitely not true. Outcomes generated by these models can be difficult to evaluate and sometimes misleading for this reason. We do not model the DCU in three-dimensions. Primarily this is to save computational resources in our simulations but most importantly a 3D model offers no further insight than a 2D model. Our 2D model effectively models the x − y coordinates of a cylindrical 3D DCU with homogeneity in the z-direction. Since we calculate the size (dimension) of this cylindrical lymph node and other pertinent parameters (such as DC flux) to be the same as its spherical counterpart, our model has the same properties as a spherical 3D model to within the accuracy of some nominal geometric prefactor. There is just one drawback to using a 2D model in this way. When modifying the DCU radius the change of volume does not scale in the same way as its 3D spherical counterpart. That being said however, in both 2D and 3D cases polynomial changes in concentration occur (as opposed to strictly linear changes if a 1D model is used) and trends in model outcomes for both a 2D model and a 3D model are qualitatively the same. We present the results from changing the ‘radius’ of the DCU in terms of fold-changes in the effective 3D volume so that it can be interpreted for both the 2D model and an equivalent 3D model. The emphasis of this manuscript is exploring important biological features of the DCU (and, by extension, the LN in general)

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

219

which play a role in the scanning of the TCR repertoire. To this end, we are more interested in trends and overall effects than a quantitatively accurate simulation (which is currently not possible with available data) and our choice to model the DCU in 2D is to capitalise on the significant computational benefits. Our argument to use a 2D model is not uncommon in this area and our justification is in line with other authors (Linderman et al., 2010). Our model consists of two cellular species, DCs and T cells. Whilst it is common for DCs in the DCU to be vastly outnumbered by T cells, we are concerned with cognate DC-T cell interaction and thus consider only a very small fraction r of T cells which are cognate for a given antigen post vaccination. Furthermore, we will consider each of the NAg cognate T cell species identically and independently and so only focus on clonotypes of a particular cognate T cell species present in the DCU (with an expected copy number of nt ). Thus from here, when we refer to T cells in the context of this model we are referring only to the clonotypes of a single T cell species, a species which is also cognate for the antigen associated with the DCs of the model. On the other hand, after a vaccination, the number of mature DCs presenting the antigen of interest is greatly increased. We therefore include the two cellular species in a hybrid model. The T cells in the model (in the DCU of a particular species) have very small copy numbers and are subject to substantial intrinsic noise. Hence we model T cells using an agent-based approach, allowing for stochastic trajectories of individual cells through the DCU. This model framework applied to DCs is computationally infeasible since the post-vaccination number of DCs is very large (on the order of 105 ). A continuum concentration described by the numerical solution of an advection-diffusion PDE is used to model the DC concentration as a function of space and time. In Appendix A, we derive the number of cognate T cell species (the cognate TCR diversity) NAg and the average number of T cell clonotypes per species Nt from T cell data in sheep as well as rates of T cell transfer between the DCU and the external lymphatic system from estimated parameters found in the literature (and presented in Table 1). The process of scanning any one species of T cell is assumed to be independent of the others and this is why we focus our model on a single species of T cell. Furthermore, only a single cell of a cognate T cell species needs to have been activated by a DC to consider that T cell species as ‘activated’. We will define an immune response to be ‘completed’ when at least one clonotype of each of the NAg cognate T cell species has been activated by interaction with a DC. In the domain  we model the DC and T cell movement assuming (in principle) both passive and active transport through the domain. This is because both DCs and T cells move through the DCU with unbiased, uncorrelated random motion as seen in two-photon experiments (Miller et al., 20 04a; 20 02), but also, on a tissue scale, move through cortical sinuses from the AL to EL (in a net direction aligned with the flow of lymph through the LN which we assume to unidirectional from top to bottom) (Sainte-Marie et al., 1982; Tilney, 1971).

parameters Dd and Vd are completely unknown at the tissue scale. In Appendix B, we derive suitable transport parameters Dd and Vd using Monte Carlo tests and matching known DC dwell times with the transport-limited removal of DCs at the medulla. We assume, however, that the dense environment of the DCU encourages diffusion-dominant transport. We set Vd = 0 for most of the simulations in the results section of this manuscript but discuss the affect of increasing active transport later in the document and so here we formally keep this functionality of the model. DCs enter the DCU at the top of the domain where the DCU interfaces with the superficial cortex; however the exact spatial distribution of inflow is unclear. In this model, we assume a uniform flux over the top of the boundary ∂ a . The extent of the afferent boundary ∂ a is described by the angle θ a which subtends it. That is, the boundary ∂ a is that part of ∂  which lies between θ = (π − θa )/2 and θ = (π + θa )/2 in standard polar coordinates (Fig. 2). Our numerical tests have indicated that the exact value of θ a has negligible consequences for the immune response and so we use a moderate value of θa = π /4 for all of our simulations (see supplementary material). The total DC influx, spread evenly over ∂ a , is a function of time G(T) and depends on the details of the vaccination or infection. DCs efficiently drain into concentrated medullary sinuses at the base of the DCU on arrival and thus an absorbing boundary condition is prescribed on ∂ e (that is, U = 0). Like ∂ a , the extent of ∂ e is described by an angle θ e which subtends it. That is, the boundary ∂ e is that part of ∂  which lies between θ = (3π − θe )/2 and θ = (3π + θe )/2 in standard polar coordinates (Fig. 2). It is likely that a LN is capable of adaptively modifying the flow of lymph through the EL (for example when the LN is inflamed) and this can be effectively modelled by increasing or decreasing the parameter θ e . On the rest of the boundary ∂ ρ = ∂  \ (∂ a ∪ ∂ e ) a no flux condition is implemented. Whilst in practice cells may flow out of the boundaries, by assuming no flux the model controls the extent of outflow using just the parameter θ e (which can be thought of as a lumped parameter describing the complicated manner by which cells leave the DCU). Modelling the cellular outflow in this particular way allows us to investigate the effect of physiological ‘bottlenecking’ in the LN (in situations such as LN shut down (Haig et al., 1999)) by varying just one parameter in a heuristically meaningful way. As we have previously done for θ a , we set a moderate baseline value for θe = π /4. Whilst θ e plays a role in controlling the removal rate of DCs, we initially set a default θe = π /4 for simulations and match the baseline removal rate β d (from Table 1) by choosing suitable transport parameters as previously stated (see Appendix B). Once suitable transport parameters are found we assume these to be constant and vary θ e to restrict cell egress from the baseline model. In this way, the actual default value of θ e is not important because choosing a different default θ e simply changes the transport parameters of the cells (which are biologically unknown at the whole tissue scale) to match known egress rates. We also set an initial (pre-vaccination) concentration of DCs to zero. The boundary and initial conditions for the DCs are summarised as follows:

2.3. Modelling DCs

U =0

We model the DC concentration by the diffusion-advection PDE

∂ U ( R, T ) ¯ ¯ U (R, T ) −VdU (R, T )), R ∈ , 0 < T < Tmax , = ∇ · ( Dd ∇ ∂T (1) where U(R, T) is the concentration of DCs inside the DCU at posi¯ is a gradient with respect to spatial coortion R ∈  and time T, ∇ dinates, Dd is the diffusion coefficient for the DCs and Vd = −Vd jˆ is the (rather crude) flow velocity of DCs towards the medulla. The

 

on

∂ e ,

(2)



∇¯ U +

VdU jˆ · nˆ = 0 Dd

∇¯ U +

 ˜ F (T ) VdU G (T ) N jˆ · nˆ = = Dd Dd D d R 0 θa

U =0

on

∂ ρ ,

(3)

on

∂ a ,

at t = 0,

ˆ is the outward facing normal on the boundary. where n

(4)

(5)

220

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

The DC influx G(T) is presented in (4) and is determined by ˜ , the total number of post-vaccination DCs that enter parameter N the DCU, and F(T), the probability density function for the delay in time between the moment of vaccination (T = 0) and DC migration into the DCU. Since DC maturation (at the site of vaccination) and migration (to the lymph node) is a multistage process with two dominant stages we model F(T) using a Gamma distribution with a shape parameter k = 2 and scale parameter θ . The probability density for the Gamma distribution is

F (T ) =





T k−1 −T exp , θ

(k )θ k

(6)

where is the Gamma function. In Neeland’s doctoral thesis ˜ are (Neeland, 2015), the scale parameter θ and total DC inflow N fit to experimental data in an ovine model for various vaccination ˜ and θ taken protocols. In this manuscript, we use the parameters N from Neeland’s thesis and presented in Table 1. These parameters correspond to fitted DC inflow data taken from a sheep receiving a standard liposomal vaccine formulation. In general, however, G(T) is a rather complex function determined by a large number of biological factors including the type and site of vaccination as well as less understood variables such as the porosity of the superficial cortex. The PDE is solved until Tmax which is determined as the time at which the vaccination is deemed to be complete. We determine this value as the time at which the total number of DCs in the DCU has reduced to a thousandth of its maximum number. 2.4. Modelling T cells A single (cognate) T cell species is described using an agentbased model. We assume T cells undergo random and directed motion in a similar fashion to DCs (albeit much more rapidly). We model individual T cell trajectories (Ri (T) for the ith T cell) through the DCU using a Langevin equation (the Fokker–Planck description of which takes the form (1) with different transport parameters)

dRi = Vt dT +



2Dt dWi ,

R∈

0 < T < Tmax

(7)

where Wi are 2D Weiner processes, Dt is the diffusion coefficient for the T cells and Vt = −Vt jˆ is the net active transport velocity of T cells towards the medulla (using similar notation as that of DCs). We track each T cell in time by updating the trajectory of each modelled T cell in the DCU after a finite time step T by taking the Ito integral of (7) from T to T + T .

R i ( T + T ) = R i ( T ) + V t T +



2Dt T ξ i ,

(8)

where each element of ξ i for each time step and each T cell are iid random numbers normally distributed with unit variance and zero mean. Just as with DCs, the transport parameters Dt and Vt are unknown but are estimated in Appendix B using known DCU dwell/exit times for T cells. T cells are supplied by HEVs into the DCU and these entry points seem to be uniformly distributed over the DCU domain. We therefore initiate T cells at random inside the DCU at a stochastic rate, per unit time, proportional to the number of circulating T cells in the lymphatic system.

(T ) = αt (Nt − n(T )),

(9)

where n is the integer valued (and piecewise constant) random variable describing the total number of modelled T cells in the DCU at time T. The parameter α t is specifically the inverse of the expected wait time to see a specific T cell enter the DCU. A more comprehensive model of the DCU may require α t , and subsequently , to depend on U (the DC concentration in the DCU) since it has been reported that DCs can stimulate an increase in transendothelial migration from HEVs into the DCU (Moussion and Girard, 2011). However, as it is unclear what this DC-dependent

control of T cell inflow should be, we ignore this affect in this model and use constant α t which is estimated based on T cell circulation time and flux τ and φ (see Appendix A and Table 1) and explore the effect of changing α t in the numerical results later. Under equilibrium conditions E (n(T )) = nt =  Nt and the constant distribution of n is given by a binomial distribution

 

Pn =

Nt n

 

 n (1 −  )Nt −n =

Nt n

αtn βtNt −n , (αt + βt )Nt

(10)

where β t is the inverse of the dwell time for T cells to persist in the DCU in a single visit (see Appendix A). We use this distribution for n to initialise the model by sampling for n and then randomly distributing that sample of n T cells in the DCU. Whilst we do this in our simulations, we acknowledge that this is somewhat unimportant because the T cells in the DCU have rapid turnover rates on the time-scale of the immune response and the number of DCs at the start of the simulation is very small. Thus, these early time points are highly inconsequential to the overall simulation and in the time it takes DCs to enter the system in any real quantities, the simulation for T cells has already reached a steady state (in an average sense). At the boundaries of the domain we have assumed no net flow of T cells, except at the boundary ∂ e where the T cells are absorbed. T cell reflection at the boundary can be easily done each time step by replacing Ri (T + T ) whenever |Ri | > R0 by (2R0 − |Ri | )Rˆ i , where Rˆ i = Ri /|Ri | unless the polar orientation of Ri lies within the interval ((3π − θe )/2, (3π + θe )/2 ) in which case it is removed from the DCU (and recirculated into the lymphatic system). Two-photon experiments have shown that T cells are able to exit through cortical sinuses spanning through the paracortical area (Grigorova et al., 2010). The cortical sinuses direct flow of T Cells from the DCU but precisely how this is done is not entirely clear. For example, some cells have been seen probing the surface for extended periods of time rather than simply entering the sinus and leaving the LN (Grigorova et al., 2009). We will model the effect of cortical sinus removal of T Cells by introducing a net active flux Vt rather than simple volumetric removal of T Cells in the DCU. The combination of active and passive transport through the DCU is controlled by total dwell time of T Cells (see Appendix B) but we explore the contributions of varying the balance of these transport modes. Determining an activation event In each time step of the model T, each T cell in the DCU is tested to see if it (and therefore the whole T cell species it belongs to) is activated by a DC. The activation rate (per unit time) for a single T cell κ is proportional to the local DC concentration U(Ri , T).

κ (Ri , T ) = U (Ri , T ).

(11)

The probability in each time step that a T cell will be scanned is given by κ (Ri , T)T where Ri is the position of the T cell at time T. We derive  for our spatially heterogeneous 2D model by adapting calculations used for the well-mixed model of Neeland (2015). In Neeland (2015), a DC-T cell interaction rate is theoretically derived for a well-mixed model of sheep LNs. This derivation is based on a Smoluchowski diffusion-limited interaction rate and is dependent on motility of both DCs and T cells (at the cellular scale found from two-photon experiments), dendrite lengths (the radius of influence of DCs in finding nearby T cells) as well as the volume of the domain (V = 1.9 × 1012 μm3 ). Here the assumption is that DCs are very efficient at interacting with nearby T cells and the limiting factor for a specific T cell finding a DC is their proximity to each other. In Neeland’s model an approximation of LN = 4π RDC D/V = 8 × 10−7 per hour is quoted as the rate

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

at which a single DC will find a particular T cell in the well-mixed domain. Here RDC is the radius of influence of DC dendrites and D is the sum of cellular scale random motility of DCs and T cells respectively. In Section 2.2, we have approximated the volume of a DCU to be a third of the volume of the domain in Neeland (2015). Thus, for a well-mixed DCU we expect a T cell-DC interaction rate of DCU = 2.4 × 10−6 per hour. We now focus on translating this rate into a 2D rate constant  given by Eq. (11). Since the rate DCU assumes the domain to be well-mixed, each T cell experiences a uniform concentration of DCs U0 (T). The rate of a particular T cell being scanned by one of nd (T) DCs in the domain  is given by κ (T ) = DCU nd (T ) = π R20 DCU U0 (T ). Thus, we have a relationship between the rate of scanning of a T cell as a function of the DC density (as required by Eq. (11)). We translate this into our spatially heterogeneous model by making the final assertion that T cells are scanned at a rate proportional to the local DC concentration (irrespective of whether it is spatially heterogeneous or not) −1 and thus we have  = π R20 DCU = 211.8 h μm2 . Occasionally, all of the DCs of the vaccination pass through the DCU without finding a clonotype of a specific T cell species. That is, the simulation makes it to Tmax without a single T cell-DC interaction. In the case where this T cell species is one of the NAg cognate species, we consider the vaccination a failure after this time. This, however, was rare with the default parameters which we were using (parameters for known successful vaccination protocols). Simulating 1.5 × 104 species of T cell experiencing a vaccination event with parameters given in Table 1 resulted in only 23 species that went unactivated by the time the DCs had entirely left the DCU. Thus the probability that a TCR of any specific type is not scanned is p ≈ 0.0015. The probability that the whole repertoire is therefore scanned (in its entirety, and given each TCR is scanned independently) is given by (1 − p)NAg ≈ 0.97. Thus,there is only a 3% chance that a vaccination of this type, drained through a single DCU, will not scan the entire cognate repertoire. In the next section we will discuss what is meant by scanning the cognate repertoire or ‘completing the immune response’ in the context of the simulations. 2.5. Determining the completion of immune response The completion of the immune response occurs when each of the NAg cognate species of T cell are activated. We assume that each of these T cell species are scanned independently. It is unclear, however, how appropriate this assumption is due to the fact that after the activation of the first T cell species the LN may undergo physiological changes before the scanning of the full TCR repertoire is achieved. Nevertheless, we interpret the time for a complete immune response under constant physiological conditions (scanning of the TCR repertoire) by simulating NAg realisations of the T cell-DC hybrid model described here. The simulated time Tˇ for the complete immune response is the time at which the last of these T cell species is activated

Tˇ = max T j , j

(12)

where Tj is the time at which the jth T cell species was simulated to have been activated, j = 1, 2, . . . , NAg . In the numerical results section of this paper we run Nsim = 150 0 0 Monte Carlo realisations of a given T cell species for any one set of parameters. After storing each of the scan times Tj we draw 10 0 0 samples of NAg = 20 with bootstrapping and present violin plots of the distribution of Tˇ (as determined by (12)) for the 10 0 0 samples. This violin plot indicates the distribution of immune response times (given an immune response). Separately, the samples which have at least one unscanned T cell species are counted and used to determine the probability of vaccination failure (that is, the probability that not

221

all cognate T cell species are activated by the DCs within a single DCU).

2.6. Practical aspects for running the simulations Since we are interested in investigating parameters which may impact domain properties we find it necessary to nondimensionalise the model described above. The details of this nondimensionalisation and other important numerical considerations such as how the T cells sample the DC distribution each time step are presented in Appendix C.

3. Results and discussion In this section we subdivide the results and discuss the sensitivity of immune response to changes in various physiological parameters. For each set of parameter sensitivity tests, each parameter is varied by various factors (given on the x-axis of each figure). The main graph shows the distribution of immune response times Tˇ using violin plots. For reference, it is important to remember that DCs are not dumped into the simulation all at once. We ˜ over time using the distribudistribute the inflow of total DCs N tion F(T) which is shown on the far left of each set of results. In all cases, an immune response cannot proceed until there are a significant number of mature DCs in the DCU so this distribution provides a soft lower bound to the expected immune response times. Above many of figures a color coded bar plot indicates, for each parameter regime, the probability that the vaccination fails. This bar graph is not included if all simulations were void of failed immune responses. Colors in each of the figures link results between the violin plots, bar graphs and DC distributions where appropriate (non-zero).

3.1. Changes in the size of the DCU In Fig. 3, it is important to note that as we change R, other parameters such as cell transport rates are being held constant. Thus, whilst the density of cells becomes smaller, we also notice that cells take longer to migrate to the medulla for removal and stay in the LN for longer. For very small DCUs under these conditions (that is, all other parameters held constant), the DCs and T cells find the EL quickly and for the entire duration of the vaccination there are very few DCs and/or T cells in the DCU at all making it rapidly unlikely (as R0 decreases) for an immune response to be mounted. Thus, when the DCU is a quarter of its estimated size a total of 0 out 10 0 0 samples reached an immune response before all of the DCs had arrived at the medulla. Increasing the size of DCUs to slightly larger than the default size for the default parameters (determined by experiments) the vaccination success rate rapidly improves to 100 percent but the actual time taken for a vaccination event to occur remains surprisingly steady (only a weak increase in time as the DCU radius becomes very large). The increased response time due to an inflammation leading to 49 times the volume of the default DCU only increases the immune response time by 10 hours and is comparable to the standard deviation in response time due purely to intrinsic stochastic fluctuations. This is due to the fact that although DCs are at a lower concentration they stay in the DCU for longer and this is the same for TCs. This numerical experiment demonstrates that the trade-off between dilution and increased dwell-time means that DCU size has a positive effect on total immune success rate (supporting the idea that inflammation leading to enlarged LNs can improve the immune response) with a minimal effect (albeit a weak general increase) on response time.

222

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

Fig. 3. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, the volume of the DCU is changed from 0.1 of its default value to 49 times its default value with very little change in response time but rapid improvement in the total probability of a successful TCR repertoire scan. Above the violin plot is a histogram of the probability that the vaccination will not elicit a full TCR repertoire scan. These bars are color coded with the violin plots since the violin plots only show the distribution of immune response times in the event of a successful immune response. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference.

Table 1 A summary of the parameter values used in our model of ovine DCUs. The table displays the mathematical symbol, the name and default value of each parameter as well as a source reference to where it has been quoted in the literature or inferred in this manuscript from the model. All figures in the results section of this manuscript are generated using the values of the parameters presented in this table, unless otherwise stated. The table is sorted by gathering related parameters. Biologically/physiologically relevant parameters regarding T cell properties, DC properties, the vaccination process and the DCU geometry are separated as well as inferred parameters which are directly used in the model. Symbol

Parameter

Value

Source

72 h 8.9 × 106 h−1 6.3 × 106 10−6 2 × 107

Thomas et al. (2012) Neeland (2015) Neeland (2015) Mirsky et al. (2011) Mirsky et al. (2011)

2 7.03 h 1.57 × 105

Neeland (2015) Neeland (2015) Neeland (2015)

1/60 h−1 211.8 h−1 μm2

Miller et al. (2004a) Neeland (2015)

π /4 π /4

5.3 × 10 μm

Section 2.2 Section 2.2 Section 2.2

1.4 × 10−2 h−1 1.4 h−1 20 32 8.2 × 105 μm2 h−1 0 μm h−1 6.9 × 107 μm2 h−1 0 μm h−1

Appendix Appendix Appendix Appendix Appendix Appendix Appendix Appendix

T cell

τ φ

T cell recirculation time T cell flux in a DCU Number of T cells in the DCU ntotal r Frequency of cognate T cells TCR diversity NTCR Vaccination k Shape parameter θ Scale parameter ˜ N Total number of DCs DC βd Rate of DC removal  T cell-DC activation rate DCU θe EL angle θa AL angle R0 DCU radius Model Parameters αt T cell input rate per T cell βt T cell removal rate per T cell Cognate TCR diversity NAg T cell Clonotypes per species Nt DC diffusion coefficient Dd DC advection Vd T cell diffusion coefficient Dt T cell advection Vt

3

A A A A B B B B

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

223

Fig. 4. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, the T cell inflow into the DCU is increased from its default value to 4 times its default value to investigate the affect of increased HEV transendothelial migration. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference.

3.2. Increasing T cell supply to the DCU It seems clear from Fig. 4 that there is a clear purpose to increasing T cell inflow during a vaccination event. Increasing the inflow by a factor of 4 substantially decreases the time for the immune response. The shortened range of the immune response time would also indicate that T cells are rapidly finding DCs after they enter the DCU. That is, the process becomes more concentrationlimited rather than transport-limited. It is important here to acknowledge that the isolated increase in T cell inflow is purely a theoretical thought experiment (much like the majority of the results in this paper) as it is difficult, in practise, to decouple this parameter from other parameters. For example, a higher T cell inflow (keeping constant the dwell time of T cells by keeping constant the transport parameters) may well increase the volume of the DCU. Since we have already seen that increasing the size of the DCU seems to have only a weak effect on the immune response time, it is clear that T cell inflow is an important factor in a successful and timely immune response. 3.3. Reduced EL lymphocyte output Reducing the egress of cells from the LN can significantly improve the immune response time. In Fig. 5, constriction of the EL and medullary sinuses by a factor of 8 seems to be all that is required to achieve optimal immune response. This optimum is the optimum in which DC and T cell egress can be ignored and the limiting factor for an immune response is the inflow of cognate T cells and mature DCs. As it is with increasing T cell inflow, decreasing cellular outflow (whilst keeping all other factors constant) would, in practise, inflame the LN (as is the case during LN shutdown). According to our 2D simulations, the combination of increased T cell inflow and reduced cellular outflow during a vaccination has a stronger effect on the immune response time than increasing the size of the LN. Our 2D simulations are equivalent to a 3D lymph node expanding in both x and y whilst not expanding in the z direction. A three-dimensional model which increases in size isotropically is likely to have a more pronounced delay in

response time but because this delay is due to decreases in density, it is expected that the overall change in equivalent volume between the two-dimensional model and a three-dimensional model is the important parameter determining this effect. Thus, we still expect that increased TC inflow and reduced cellular outflow are the most important factors determining the response time. 3.4. The effect of macro-transport mechanism It is surprising that, for a process which is transport-limited, the mechanism of transport has such a weak (if not zero) effect on the immune response. In Fig. 6 we can see that increasing the Péclet number from 0 to 1 results in a negligible effect on the immune response. Roughly speaking the Péclet number is the ratio of diffusive transport to active transport. A more detailed description of this is given in Appendix B. Briefly, it increases the amount of motion due to active transport experienced by the cells. However, it is important to note here that by increasing advection (and therefore Péclet number) we have ensured that the overall dwell time remains constant (see Appendix B). We have already seen that decreasing the dwell time has an adverse effect on the immune response. Here we show that, for all its trafficking complexity, the precise manner by which cells go from AL to EL in the DCU is not important to the objective of scanning T cells. 3.5. Variations in DC input and vaccination protocols The size and duration of an infection/vaccination event and its impact on the immune response is an important property to investigate. Vaccine technology requires a precise dose be administered in a particular way. New technologies are making it possible to increase mature DC numbers in the LN and modulate the time scale over which they drain into the LN. We investigate effect of increasing the time scale over which the DCs enter the LN θ and the total number of mature DCs as a result of the whole vaccination event ˜ in Figs. 7 and 8 respectively. N The duration of the DC inflow does not seem to have a significant effect on the immune response time when all other param-

224

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

Fig. 5. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, the exit angle θ e is reduced from its default value of π /4 by factors 1/2, 1/4, 1/8 and 1/16 to investigate the reduction in cell egress commonly observed during an infection. Above the violin plot is a histogram of the probability that the vaccination will not elicit an immune response. These bars are color coded with the violin plots since the violin plots only show the distribution of immune response times in the event of a successful immune response and is only non-zero for the default parameters. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference.

Fig. 6. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, we increase the Péclet number (the dominance of active transport of cells over passive transport of cells) from 0 (in the default case) to 1 (where advection and diffusion are roughly equally responsible for the trafficking of the cells through the DCU). Above the violin plot is a histogram of the probability that the vaccination will not elicit an immune response. These bars are color coded with the violin plots since the violin plots only show the distribution of immune response times in the event of a successful immune response. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference.

eters are kept constant. In general, there does seem to be some increase in the response time. This is not to be unexpected since the DCs are just taking longer on average to enter the DCU. Interestingly, when the delay between vaccination and DC inflow into the LN is particularly short or long the immune response seems to be slightly less reliable than the parameters in between. Of course, these differences are miniscule and not statistically significant for the set of parameters we are using and more investigation should be made into this. It is theoretically possible however, that since cognate T cells in the DCU are highly stochastic in number and there are long periods without any cognate T cells of a particular species, it may be possible that if post-vaccination DCs drain into the DCU over a short interval of time and are removed too quickly,

the DCs may be all cleared before any T cell of a particular species visits the draining DCU. On the other hand, if the DCs drain into the DCU over a long period (without changing the total number of DCs) the DCU may be too dilute of DCs to mount an effective (timely) immune response. The impact of number of DCs is not very surprising. Increasing the number of DCs (and keeping all other parameters constant) not only speeds the immune response up but also rapidly increases the likelihood of a complete TCR repertoire scan. In practise however, a vaccine that causes the maturation of a lot of DCs may also make the patient ill and optimal vaccine doses should also include costs for likelihood of adverse reactions to the vaccine (not included in this model). This model predicts, however,

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

225

Fig. 7. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, we vary the parameter θ from 0.5 to 2.5 of its default value. Above the violin plot is a histogram of the probability that the vaccination will not elicit an immune response. These bars are color coded with the violin plots since the violin plots only show the distribution of immune response times in the event of a successful immune response. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference. Since this entry time distribution for DCs is determined exclusively by the parameter θ the distributions are also color coded to the violin plots for which they are relevant.

Fig. 8. Violin plot of immune response times Tˇ of a vaccination using parameters given in Table 1. On the horizontal axis, the total number of mature DCs which enter the DCU as a result of the vaccination is varied from 0.2 to 1.1 of default value. Above the violin plot is a histogram of the probability that the vaccination will not elicit an immune response. These bars are color coded with the violin plots since the violin plots only show the distribution of immune response times in the event of a successful immune response. Importantly, at 0.2 of the default mature DCs none of our 10 0 0 sample simulations resulted in a successful immune response which is why there is no violin plot for this value. To the left of the violin plot, a distribution for the delay time between vaccination (at time t = 0) and DC entry into the DCU is presented for reference.

226

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

˜ as inferred from experiments is the DC copy number that that N is required to push the likelihood of vaccination failure to 3% (in a single DCU). Since many of the parameters in this model have had to be inferred from literature, this is a fairly reliable/realistic result to obtain purely from simulation given the complexity of the immune response. This model could be used as an adjunct tool for optimising vaccine protocols with more precise parameters. This is made especially true since we have found that the mode of macroscopic cellular transport within the DCU environment (which is quite complex and probably the most unknown characteristic of the model) plays little to no role in the immune response. Throughout this manuscript we have analysed the factors which determine the immune response as a result of mature DCs entering a single DCU to find cognate naive TCs. This may naturally raise the question about the immune response that occurs as a result of DCs entering more than one DCU (for example, within the same draining LN or in another draining LN). Multiple DCUs have not been considered here because the variability in DCU conditions simply add to the variability in immune response. For example, consider ˜ DCs (see Table 1) which pass through a single DCU or evenly N distributed between two identical DCUs. It can be shown from Fig. 8 that if all DCs enter a single DCU then the probability that the DCU scans the full repertoire is ∼ 0.97. Alternatively, the probability that the full repertoire is scanned within a single DCU if only half of the DCs pass through it (and the other half through another DCU) is ∼ 0.5. This gives a probability of failing to scan any particular TCR in any one of the two DCUs of p = 1 − 0.51/NAg ∼ 0.034 and a probability of scanning at least one of each TCR between the two DCUs of (1 − p2 )NAg ∼ 0.97. The same as if they had both gone through the same DCU. This is because, assuming that all DCs and TCs act independently, it doesn’t matter which DCU a particular DC is drained into, it experiences the same TC environment (on an average level). Thus, all of the immune response results for a single DCU can be considered identical to those results expected if the same number of DCs were distributed among multiple DCUs. 4. Conclusion In this manuscript we present a tissue scale spatially resolved mathematical model of the LN (focusing on the DCU specifically) during a vaccination. We use a novel hybrid mathematical model which treats the large number of DCs which arrive postvaccination/post-infection as a continuum whilst naive cognate T cells of any given species is treated in a stochastic manner (one cell at a time). We base all of our parameters for this model on values from the literature or fit to data we have obtained in our lab and we focus on a sheep model. Our investigation is focused on a discussion of which physiological factors are the most important to a successful (and timely) immune response. Interestingly, the most surprising results that we found were the robustness of immune response times to many of the physiological factors and the fact that other factors seem to have little effect on the immune response in general. Cells are trafficked through the LN with elaborate micro-infrastructure. We found that, as far as efficient immune response and scanning the TCR repertoire is concerned, it does not matter if the tissue scale trafficking of the cells is dominated by active transport or passive transport but rather just the total time the cells reside in the DCU. This is surprising given the known transport-limited interaction between DCs and TCs. Keeping the transport parameters constant and increasing the radius of the DCU means that cells will reside for longer in the DCU but also that DCs will become more dilute. The overall result is a slight increase in the immune response time and improvement in overall immune success. The most important factors for a successful (and timely) immune response as found by our simulations unsurprisingly include

increased T cell supply rate, reduction in cellular egress and increased mature DC copy numbers. The later of these three factors often cannot be favourably changed without adverse effects for a patient. Since it has been shown that DCs may stimulate higher T cell inflow, our results support efforts to improve this mechanism. Our model gives reasonable results using experimentally derived parameters. It is unclear whether a more elaborate treatment of cellular flow and trafficking which correctly treats the boundary interchange of cells would offer any more insight than what has been achieved in this work since our rudimentary investigation of transport mechanisms seem to suggest no role for the nature of the flow. We believe that a model of this type holds promise as an adjunctive tool for optimising vaccinations in the future. Appendix A. Derivation of ovine T cell trafficking parameters In Appendix A, we focus on estimating important model parameters from other known parameters specific to sheep. Any given parameters in this appendix are known from the literature or from our lab experiments (see Table 1 for references). We begin with estimates of T cell types and copy numbers in sheep. The TCR diversity, the number of total T cell species, in sheep is estimated to be NT CR = 2 × 107 . We could easily find the expected number of clonotypes per species Nt if we have the total number of T cells in the whole sheep Ntotal (since Nt = Ntotal /NT CR ) but this number is not usually quoted in the literature. Instead, we infer Ntotal based on trafficking estimates and observations. It is known that T cells take an average 24 h to circulate through the entire lymphatic system (and revisit the same LN). That is, we will expect every T cell to visit our DCU about once every τ = 72 h (recalling our assumption of about 3 DCUs in a LN). We say that

τ = τDCU + τSY S =

1

βt

+

1

αt

,

(A.1)

where τ DCU and τ SYS are the expected durations of time a particular T cell spends in the DCU and ‘the rest’ of the lymphatic system respectively in that 72 h period. We will find it more convenient to use the rates (rather than durations) β t and α t at which a single T cell leaves and returns to (resp.) the DCU. On the other hand, in homeostatic conditions a DCU does not change in volume or expected number of cells and so conservation of cellular copy numbers requires that the flux φ (total number of cells per unit time) entering the DCU equal to the flux leaving the DCU. This gives us the conservation equation

φ = αt (Ntotal − ntotal ) = βt ntotal .

(A.2)

Neeland measured this flux of T cells in her thesis by measuring the number of T cells leaving the EL (of the whole LN) per unit time and found φ = (2.67 × 107 )/3 = 8.9 × 106 h−1 . Neeland also measured a total of 1.9 × 107 T cells in a LN. Assuming these T cells are found predominantly in our three DCUs, this means that the total number of T cells in a DCU is ntotal ∼ 6.3 × 106 . This allows us to calculate our DCU transit rates α t and β t . Using (A.2) we have that

βt =

φ ntotal

−1

= 1.41h

.

(A.3)

Rearranging (A.1) we have that

αt =

βt −1 = 1.40 × 10−2 h . τ βt − 1

(A.4)

Using the other equality in (A.2), (A.3) and (A.1) we can now estimate the total number of T cells in the whole lymphatic system, as expected it is just the flux multiplied by the circulation time (the time taken to count every T cell once at a point)

Ntotal =

φ φ φ + ntotal = + = τ φ = 6.41 × 108 . αt αt βt

(A.5)

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

Furthermore, we can now calculate the expected number of T cell clonotypes in each species

Nt =

Ntotal = 32.0. NT CR

(A.6)

For completeness we also define here  as the ratio between the scale of the DCU (in our model) compared to that of the whole lymphatic system.

=

ntotal = Ntotal

αt ∼ 0.01. αt + βt

(A.7)

We will be assuming that the DCU is an unbiased sample of the T cells in the whole lymphatic system and so scaling from the whole organism to our DCU just requires multiplication by  . We now shift our attention to cognate T cells. It has been observed in sheep that approximately a fraction r = 1 × 10−6 of all naive T cells will be cognate for a given antigen. Given that there are NT CR = 2 × 107 species of naive T cell, we expect

NAg = rNT CR = 20

(A.8)

of them to be cognate for the antigen being carried by the DCs. Using these parameters it is simple to calculate others. For example, nt , the expected number of T cells of a given species in the DCU at any moment in time is simply given using the ratio  ; nt =  Nt = 0.32. Importantly this means that for about two thirds of the time a given T cell species will have no representative clonotypes inside a given DCU. One would imagine then that if DCs pass through the DCU too quickly, they may not ever have the chance to encounter a T cell. It is important therefore for our model to be stochastic in the number of T cells present in the DCU. Another parameter that can be easily derived here is the total number of T cells in the DCU of any species cognate for the antigen. This is given using the ratio r; rntotal , or calculated indirectly using r; NAg nt . In this way, scaling between DCU and whole lymphatic system involves the ratio  and scaling between cognate T cells as a fraction of all T cells uses involves the ratio r. It is also important to note that these parameters are only rough parameters in order to get realistic results. It is accepted in the literature that these values have distributions associated with them and that these distributions may be substantially skewed (Asquith et al., 2015). As such, we uses these parameters only to two significant figures (see Table 1). Appendix B. Derivation of DC and T cell transport parameters Whilst there is a fair body of literature that focuses on the transport of DCs and T cells at the cellular scale (mostly using twophoton experimental data), tissue scale cellular transport parameters Dd , Vd , Dt and Vt are completely unknown. We will not attempt to presume calculating the tissue scale transport parameters for DCs and T cells from cellular processes but instead take a topdown approach. It is stressed for the reader that these transport parameters at the tissue scale are not the same as those that one might find at the microscale using, for example, two-photon experiments. The transport parameters, which we will fit to known macroscopic cellular trafficking rates, are homogenised multiscale parameters heuristically taking into consideration all complex interactions inside the DCU. We use known ovine DC removal rates βd = 1/τd (see Table 1) and T cell removal rates βt = 1/τt (derived from T cell flux data in Appendix A). These removal rates define expected dwell times for both DCs and T cells τ d and τ t (the expected times these cells persist in the DCU). Since the dwell times are controlled by the transport parameters (that is, they are removed only when they reach the EL) we can use these dwell times to infer transport parameters (Dd , Vd , Dt and Vt ). Since there is too many degrees of freedom in the parameters Dd , Vd , Dt and Vt to

227

get a unique solution for a given τ d and τ t , we set Vd and Vt and calculate Dd and Dt . As we are unsure to what extent active transport plays in trafficking cells through the DCU we initially set Vd = Vt = 0 and later investigate the qualitative affect on the immune response in the presence of some active transport (but keeping dwell times constant). The calculation of dwell times for Brownian molecules inside circular (spherical) domains with a designated target/exit is something that has been studied extensively in applied mathematics (and in particular mathematical biology). Such problems are notoriously difficult to obtain closed-form solutions for and there exist only closed-form solutions for the dwell time with specific conditions such as in infinite domains and spherical targets (Smoluchowski, 1918) and in rounded cavities with small exits (Taflia and Holcman, 2007). We shall use Monte Carlo simulations to tabulate dwell times. Consider a cell (either a T cell or a DC) with a diffusion constant D and velocity −V ˆj inside a circle of radius R0 and a constant EL aperture angle θe = π /4 describing the domain  outlined in Section 2.3 (with reflective boundaries outside of ∂ e ). We denote the expected time that this cell spends inside the circle τ exit ≡ τ exit (D, V, R0 , ρ), where ρ is its initial position inside the cell. We wish to find D given R0 , V, and τ dwell (where τ dwell is the known dwell time for that cellular species: τ d for DCs and τ t for T cells) such that

Texit (D, R0, V ) =





P0 (ρ )τexit (D, V, R0 , ρ )dV = τdwell ,

(B.1)

where Texit is the expected exit time for a cell which has a starting position which is randomly distributed. Here, P0 (ρ) is the probability density of the starting position of the cell (uniformly distributed over ∂ a for DCs and uniformly distributed over  for T cells) and the integral is taken over all starting positions ρ. Non-dimensionalisation of (B.1), scaling spatial dimensions by R0 and time dimensions by R20 /D gives

T¯exit (V¯ ) =



ω

P¯0 (ρ¯ )τ¯exit (1, V¯ , 1, ρ¯ )dV¯ = τ¯dwell ,

(B.2)

where the bar notation indicates the dimensionless associated variable/parameter (noting that both D and R0 have a dimensionless unit value) and ω is the unit circle. Here the non-dimensionalised velocity V¯ is the Péclet number associated with the length-scale of the DCU. To solve (B.2) (and by extension (B.1)) independently for DCs and T cells, we run Monte Carlo simulations to obtain T¯exit (V¯ ) in each case. We recall that T¯exit (V¯ ) is the expected time for a Brownian trajectory, with unit diffusion constant and given drift velocity −V¯ ˆj, starting at a random location with distribution P¯0 (which is uniformly distributed on ωa for DCs and within ω for T cells) to reach ωe . We therefore calculate this expected time by simulating these cells M times (resampling the starting position each time) and finding the average exit time - the average time it takes for the trajectory to get to ωe . To simulate the trajectories we use a generalised form of (8)

r(t + t ) = r(t ) +

√ ¯ t. 2t ξ + V

(B.3)

The error associated with calculating T¯exit (V¯ ) in this way will rely on having a large M and a small t (compared with the time-scale of the transport). In all of our simulations we chose 2 / (2 )), where  −3 M = 104 and t = min(t /V¯ ,  t = 5 × 10 . The t time-step t is chosen such that the expected Brownian displacement and expected displacement due to drift each time-step is bounded by  t and the accuracy of the Monte Carlo simulations is therefore controlled in the case of both diffusion-dominated and advection-dominated transport.

228

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

2.5

2

Dendritic cells (DCs) 1.5

T cells (TCs) 1

0.5

0 0

1

2

3

4

5

6

7

8

9

10

Fig. B.9. The dimensionless dwell time T¯ using Monte Carlo simulation for both DCs (entering at the AL boundary) and T cells (entering uniformly on the domain) versus increasing Péclet number V¯ and fit to the polynomial fraction of the form (B.4).

Table B.2 Polynomial fraction model fits to Monte Carlo simulations of T¯exit (see Eq. (B.4)). The error margins are given as a 95% confidence interval. Model parameter

DC value

T cell value

p1 p2 q1 q2

2.03 ± 0.04 6.02 ± 0.8 2.55 ± 0.15 2.62 ± 0.34

1.086 ± 0.03 4.25 ± 0.54 2.34 ± 0.22 2.436 ± 0.27

The dimensionless exit time T¯exit depends on the Péclet number V¯ but also the distribution of initial positions P¯0 . Fig. B.9 displays the dimensionless exit times T¯exit (V¯ ) with θe = π /4 and θa = π /4 calculated using Monte Carlo simulations. In order to plot T¯exit as a function of V¯ we took averages of the M Monte Carlo simulations for each V¯ ∈ [0, 10] on intervals of V¯ = 0.1. We then fit these averages to a polynomial fraction of the form



T¯exit V¯ =

p1V¯ + p2 . V¯ 2 + q1V¯ + q2

(B.4)

We found a model of this form to very closely match the data with a RMSE of 5 × 10−3 for DC simulations (simulations that start on ωa ) and 4 × 10−3 for T cell simulations (simulations that start on uniformly at random inside ω). Table B.2 summarises the model fits for both T cells and DCs. We use Eq. (B.2) to calculate D given a chosen value of V. To do this, we either assume a given Péclet number V¯ and redimensionalise Eq. (B.2) to find D



D=

T¯exit V¯ R20

τdwell

.

(B.5)

Alternatively, if the velocity is known, we can find D by solving the implicit equation

D=

R20

τdwell

T¯exit

V R  0

D

.

(B.6)

In this case, we can non-dimensionalise D and V with respect to time-scale τ dwell and length-scale R0 . In this non-dimensionalisation (denoted by a hat) Dˆ = Dτdwell /R20 , Vˆ = V τdwell /R0 and

 

Dˆ = T¯exit

Vˆ . Dˆ

(B.7)

We can determine D = Vˆ R20 /(τdwell xˆ) (where xˆ = Vˆ /Dˆ ) using (B.7) and fixed point iterations. Starting at xˆ0 = Vˆ /T¯exit (0 ) we iterate using the following fixed point relation xˆn+1 = f (Vˆ , xˆn ) = Vˆ /T¯exit xˆn . Using the model (B.4) and the parameters in the Table B.2, we can see that f is continuous and differentiable for xˆ ∈ (0, ∞ ). Furthermore, f is positive and f is also positive since for both curves p2 q1 > p1 q1 and so is f since p1 p2 q1 < p21 q2 + p22 . Thus, if the limit exists, the derivative is bounded by | f (x )| < limx→∞ f (x ) = Vˆ /p1 . From the existence and uniqueness theorem for fixed point iterations, we expect convergence of the fixed point iterations to xˆ = Vˆ /Dˆ which is the unique implicit solution to (B.7) for 0 < Vˆ < p1 but find that convergence is achieved with our seed value for an even wider range: 0 < Vˆ < 2.066 for DCs and 0 < Vˆ < 1.279 for T cells. Values of Vˆ outside this critical range are not physical since for velocities above these intervals there is no solution to (B.7). The physical basis for this is that Vˆ is the ratio between the speed of the advection and the speed required to move a distance of the radius of the circle in the dwell time. Thus, we can interpret the critical dimensionless velocities Vˆ here (2.066 for DCs and 1.279 for T cells) as the dimensionless velocity required to get from the initial positions to the EL in the dwell time (which is 1 when non-dimensionalised). Since the DCs start approximately 2 radii from the EL and the T cells on average start 1 radii from the EL these critical velocities are approximately 2 and 1 respectively. The discrepancy from 2 and 1 can be explained due to increased time that is required for some cells to ‘slide’ down the bottom of the DCU either side of ωe . The relationship Dˆ (Vˆ ), given implicitly by (B.7) and θa = θe = π /4 for both DCs and T cells is presented in Fig. B.10. These lines can be thought of as lines of constant dwell time in D − V parameter space given DC and T cell starting locations and DCU radius. Our choice of D in our model is dictated by this curve irrespective of changes in parameters such as R0 in simulations. This is because we assume D and V to be ‘fundamental’ constants which are unchanged if the LN size is changed. The dwell time that is used to generate these curves is estimated for a standard LN and this is why we use a constant R0 in the calculation of D and V. For the bulk of the manuscript the advection is considered small and thus we choose V = V¯ = 0 and D is given simply by substituting V¯ = 0 into Eq. (B.5). This is the value that, for the given R0 and β = 1/τdwell , is quoted in Table 1 for both DCs and T cells.

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

229

2.5

2

Dendritic cells (DCs) T cells (TCs)

1.5

1

0.5

0 0

0.5

1

1.5

2

Fig. B.10. Normalised diffusion Dˆ versus normalised advection velocity Vˆ which solves the implicit relation (B.7) for both T cells and DCs. That is, the relationship between diffusion and velocity such that the dwell time remains constant (and normalised to 1 for a DCU of unit radius).

Appendix C. Simulations and non-dimensionalisation Whilst we report all of our results using dimensionalised quantities, we run our simulations on the same domain and mesh and thus we non-dimensionalise. We non-dimensionalise all of our simulations by scaling space with respect to R0 and time with respect to R20 /Dd (the characteristic time for DCs to traverse distances on the scale of the DCU). The non-dimensional variables for position, time and concentration are indicated using lower case variables associated with their dimensional counterparts (including all symbols with markings and indices as well as Greek letters). For example, non-dimensional immune response time would become tˇ rather than the dimensional Tˇ . We have therefore

R = R0 r;

T =

R20 t; Dd

U=

1 u. R20

(C.1)

Where dimensional parameters are already denoted using lower case characters (including Greek letters), we have placed a bar above them to indicate the non-dimensional analogues (such as α¯ t instead of α t ). The only remaining exception to this notation is the non-dimensionalised diffusion constant for T cells which is denoted γ (and is the ratio between T cell and DC diffusion constants). The non-dimensionalised domain of the DCU ω is a unit circle. DC concentration on ω is solved using the nondimensionalised Eq. (1) and relevant boundary and initial conditions

∂ u ( r, t ) = ∇ · (∇ u(r, t ) − vd u(r, t )), r ∈ ω, 0 < t < tmax , (C.2) ∂t where vd = Vd R0 /Dd and the associated boundary conditions are

u=0

on

∂ωe ,



∇ u + vd ujˆ · nˆ = 0



∇ u + vd ujˆ · nˆ = g(t )

u=0

(C.3)

on

∂ωρ , on

∂ωa ,

at t = 0 ˜ F (R2 t/Dd )/(Dd θa ) R20 N 0

(C.4)

(C.5)

(C.6)

where g(t ) = and vd = Vd R0 /Dd is the DC Péclet number. The non-dimensionalised equation for DC concentration (C.2) is solved using the inbuilt finite element solver in MATLAB’s PDE toolbox, subject to the boundary conditions and initial conditions in Eqs (C.3)–(C.6).

Table C.3 A summary of the non-dimensionalised parameter values used in the (default) model simulations. It is important to note that R0 is crucial to the non-dimensionalisation (as it is involved in the characteristic time-scale and length-scale) and these numbers change when exploring the parameter space. Non-dimensionalised parameter

Value

Definition

γ

63.8 (0,0) (0,0) 0.37 2.0 × 10−4

D t /D d Vd R0 /Dd Vt R0 /D d αt R20 /Dd R20 /Dd

vd vt α¯

λ

T cells are updated using the following non-dimensional form of (8)

ri (t + t ) = ri (t ) + vt t +



2γ t ξ i ,

(C.7)

where vt = Vt R0 /Dd . T cells are introduced into the DCU at a nondimensional rate of

σ (t ) = α¯ t (Nt − n(T )),

(C.8)

being added with a probability of σ t each time step, and activated by DCs at a non-dimensional rate of

κ¯ (ri , t ) = λu(ri , t ),

(C.9)

activating with a probability of κ¯ t each time step. In order that the discrete solution to the Langevin equation (Eq. (C.7)) and the treatment of boundaries as described in main text is accurate, we reduce t so that the RMS of |ri (t + t ) − ri (t )| is bound by δ r = 0.01 (significantly smaller than the features of ω). Ensuring the accuracy of T cell trajectories requires a very small time-step. By a substantial margin the most time consuming step in the simulations is finding (or interpolating) u(ri , t) at each time step for use in calculating the probability of an activation event. Prior to running the T cell simulations, we tabulate the DC concentration u(r, t) (from Eq. (C.2)) for 10 0 0 time points on the interval (0, tmax ) and interpolated onto a 2D Cartesian grid (xm , ym ) where xm = −1 + mδ x and ym = −1 + mδ y where δx = δy = 0.01 and m = 0, 1, . . . , 2/δx . This grid is sufficiently fine that u(ri , t) is estimated at each time-step with a linear interpolation between time points of the tabulated DC concentration at the nearest grid point (xm , ym ) to ri . Estimating u(ri , t) in this way (rather than interpolating using the triangular mesh) allows for immediate evaluation after quickly rounding ri to the nearest grid point (ultimately speeding simulations up by a factor greater than 100).

230

D.L. Maderazo et al. / Journal of Theoretical Biology 454 (2018) 215–230

In Table C.3 we present a list of important non-dimensionalised parameters whose value (corresponding to the dimensional counterparts in Table 1) are also displayed. Supplementary material Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.jtbi.2018.06.008. References von Andrian, U.H., Mempel, T.R., 2003. Homing and cellular traffic in lymph nodes. Nat. Rev. Immunol. 3 (11), 867. Asquith, R., Laydon, D., Bangham, C., 2015. Estimating t-cell repertoire diversity: limitations of classical estimators and a new approach. Bains, I., Antia, R., Callard, R., Yates, A.J., 2009. Quantifying the development of the peripheral naive cd4+ t-cell pool in humans. Blood 113 (22), 5480–5487. Bajaria, S.H., Webb, G., Cloyd, M., Kirschner, D.E., 2002. Dynamics of naive and memory cd4+ t lymphocytes in hiv-1 disease progression. Bauer, A.L., Beauchemin, C.A., Perelson, A.S., 2009. Agent-based modeling of host— pathogen systems: the successes and challenges. Inf. Sci. 179 (10), 1379–1389. Beltman, J.B., Marée, A.F., De Boer, R.J., 2007a. Spatial modelling of brief and long interactions between t cells and dendritic cells. Immunol. Cell Biol. 85 (4), 306. Beltman, J.B., Marée, A.F., Lynch, J.N., Miller, M.J., de Boer, R.J., 2007b. Lymph node topology dictates t cell migration behavior. J. Exp. Med. 204 (4), 771–780. Bogle, G., Dunbar, P., 2010a. Agent-based simulation of t-cell activation and proliferation within a lymph node. Immunol. Cell Biol. 88 (2), 172–179. Bogle, G., Dunbar, P., 2010b. T cell responses in lymph nodes. Wiley Interdiscip. Rev. 2 (1), 107–116. Bousso, P., 2008. T-cell activation by dendritic cells in the lymph node: lessons from the movies. Nat. Rev. Immunol. 8 (9), 675. Bousso, P., Robey, E., 2003. Dynamics of cd8+ t cell priming by dendritic cells in intact lymph nodes. Nat. Immunol. 4 (6), 579. Casrouge, A., Beaudoing, E., Dalle, S., Pannetier, C., Kanellopoulos, J., Kourilsky, P., 20 0 0. Size estimate of the αβ tcr repertoire of naive mouse splenocytes. J. Immunol. 164 (11), 5782–5787. Celli, S., Day, M., Müller, A.J., Molina-Paris, C., Lythe, G., Bousso, P., 2012. How many dendritic cells are required to initiate a t-cell response? Blood 120 (19), 3945–3948. Cyster, J.G., Schwab, S.R., 2012. Sphingosine-1-phosphate and lymphocyte egress from lymphoid organs. Annu. Rev. Immunol. 30, 69–94. Förster, R., Braun, A., Worbs, T., 2012. Lymph node homing of t cells and dendritic cells via afferent lymphatics. Trends Immunol. 33 (6), 271–280. Gong, C., Mattila, J.T., Miller, M., Flynn, J.L., Linderman, J.J., Kirschner, D., 2013. Predicting lymph node output efficiency using systems biology. J. Theor. Biol. 335, 169–184. Grigorova, I.L., Panteleev, M., Cyster, J.G., 2010. Lymph node cortical sinus organization and relationship to lymphocyte egress dynamics and antigen exposure. Proc. Nat. Acad. Sci. 107 (47), 20447–20452. Grigorova, I.L., Schwab, S.R., Phan, T.G., Pham, T.H., Okada, T., Cyster, J.G., 2009. Cortical sinus probing, s1p 1-dependent entry and flow-based capture of egressing t cells. Nat. Immunol. 10 (1), 58. Haig, D., Hopkins, J., Miller, H., 1999. Local immune responses in afferent and efferent lymph. Immunology 96 (2), 155. Hopkins, J., McConnell, I., Pearson, J., 1981. Lymphocyte traffic through antigen-stimulated lymph nodes. ii. role of prostaglandin e2 as a mediator of cell shutdown.. Immunology 42 (2), 225. Kindt, T.J., Goldsby, R.A., Osborne, B.A., Kuby, J., 2007. Kuby immunology. Macmillan. Lee, M., Mandl, J.N., Germain, R.N., Yates, A.J., 2012. The race for the prize: T-cell trafficking strategies for optimal surveillance. Blood 120 (7), 1432–1438. Linderman, J.J., Riggs, T., Pande, M., Miller, M., Marino, S., Kirschner, D.E., 2010. Characterizing the dynamics of cd4+ t cell priming within a lymph node. J. Immunol. 184 (6), 2873–2885.

Marino, S., Pawar, S., Fuller, C.L., Reinhart, T.A., Flynn, J.L., Kirschner, D.E., 2004. Dendritic cell trafficking and antigen presentation in the human immune response to mycobacterium tuberculosis. J. Immunol. 173 (1), 494–506. Meyer-Hermann, M.E., Maini, P.K., 2005. Interpreting two-photon imaging data of lymphocyte motility. Phys. Rev. E 71 (6), 061912. Miao, H., Hollenbaugh, J.A., Zand, M.S., Holden-Wiltse, J., Mosmann, T.R., Perelson, A.S., Wu, H., Topham, D.J., 2010. Quantifying the early immune response and adaptive immune response kinetics in mice infected with influenza a virus. J. Virol. 84 (13), 6687–6698. Miller, M.J., Hejazi, A.S., Wei, S.H., Cahalan, M.D., Parker, I., 2004a. T cell repertoire scanning is promoted by dynamic dendritic cell behavior and random t cell motility in the lymph node. Proc. Natl. Acad. Sci. U.S.A. 101 (4), 998–1003. Miller, M.J., Safrina, O., Parker, I., Cahalan, M.D., 2004b. Imaging the single cell dynamics of cd4+ t cell activation by dendritic cells in lymph nodes. J. Exp. Med. 200 (7), 847–856. Miller, M.J., Wei, S.H., Parker, I., Cahalan, M.D., 2002. Two-photon imaging of lymphocyte motility and antigen response in intact lymph node. Science 296 (5574), 1869–1873. Mirsky, H.P., Miller, M.J., Linderman, J.J., Kirschner, D.E., 2011. Systems biology approaches for understanding cellular mechanisms of immunity in lymph nodes during infection. J. Theor. Biol. 287, 160–170. Moussion, C., Girard, J.-P., 2011. Dendritic cells control lymphocyte entry to lymph nodes through high endothelial venules. Nature 479 (7374), 542. Murphy, K., Travers, P., Walport, M., 2007. Principles of innate and adaptive immunity. Janeway’s Immunobiol. 1–38. Naylor, K., Li, G., Vallejo, A.N., Lee, W.-W., Koetz, K., Bryl, E., Witkowski, J., Fulbright, J., Weyand, C.M., Goronzy, J.J., 2005. The influence of age on t cell generation and tcr diversity. J. Immunol. 174 (11), 7446–7452. Neeland, M., de Veer, M., Scheerlinck, J.-P., 2017. The role of antigen presentation and innate immunity during immune induction by particulate antigens. In: Micro and Nanotechnology in Vaccine Development. Elsevier, pp. 83–98. Neeland, M.R., 2015. Characterising the Immune Response to Liposomal Adjuvant Formulations: A Systems Biology Approach. Ph.D. thesis. Monash University. Neeland, M.R., Elhay, M.J., Nathanielsz, J., Meeusen, E.N., de Veer, M.J., 2014. Incorporation of cpg into a liposomal vaccine formulation increases the maturation of antigen-loaded dendritic cells and monocytes to improve local and systemic immunity. J. Immunol. 192 (8), 3666–3675. Ohtani, O., Ohtani, Y., 2008. Structure and function of rat lymph nodes. Arch. Histol. Cytol. 71 (2), 69–76. Randolph, G.J., Angeli, V., Swartz, M.A., 2005. Dendritic-cell trafficking to lymph nodes through lymphatic vessels. Nat. Rev. Immunol. 5 (8), 617. Riggs, T., Walts, A., Perry, N., Bickle, L., Lynch, J.N., Myers, A., Flynn, J., Linderman, J.J., Miller, M.J., Kirschner, D.E., 2008. A comparison of random vs. chemotaxis-driven contacts of t cells with dendritic cells during repertoire scanning. J. Theor. Biol. 250 (4), 732–751. Sage, P.T., Carman, C.V., 2009. Settings and mechanisms for trans-cellular diapedesis. Front. Biosci. 14, 5066. Sainte-Marie, G., Peng, F.-S., Belisle, C., 1982. Overall architecture and pattern of lymph flow in the rat lymph node. Dev. Dyn. 164 (4), 275–309. Smoluchowski, M.v., 1918. Versuch einer mathematischen theorie der koagulationskinetik kolloider lösungen. Zeitschrift für physikalische Chemie 92 (1), 129–168. Taflia, A., Holcman, D., 2007. Dwell time of a brownian molecule in a microdomain with traps and a small hole on the boundary. J. Chem. Phys. 126 (23), 234107. Thomas, N., Matejovicova, L., Srikusalanukul, W., Shawe-Taylor, J., Chain, B., 2012. Directional migration of recirculating lymphocytes through lymph nodes via random walks. PLoS ONE 7 (9), e45262. Tilney, N.L., 1971. Patterns of lymphatic drainage in the adult laboratory rat.. J. Anat. 109 (Pt 3), 369. Tomura, M., Hata, A., Matsuoka, S., Shand, F.H., Nakanishi, Y., Ikebuchi, R., Ueha, S., Tsutsui, H., Inaba, K., Matsushima, K., et al., 2014. Tracking and quantification of dendritic cell migration and antigen trafficking between the skin and lymph nodes. Sci. Rep. 4, 6030. Willard-Mack, C.L., 2006. Normal structure, function, and histology of lymph nodes. Toxicol. Pathol. 34 (5), 409–424.