Systematic Screening of Chemokines to Identify Candidates to Model and Create Ectopic Lymph Node Structures for Cancer Immunotherapy

The induction of ectopic lymph node structures (ELNs) holds great promise to augment immunotherapy against multiple cancers including metastatic melanoma, in which ELN formation has been associated with a unique immune-related gene expression signature composed of distinct chemokines. To investigate the therapeutic potential of ELNs induction, preclinical models of ELNs are needed for interrogation of these chemokines. Computational models provide a non-invasive, cost-effective method to investigate leukocyte trafficking in the tumor microenvironment, but parameterizing such models is difficult due to differing assay conditions and contexts among the literature. To better achieve this, we systematically performed microchemotaxis assays on purified immune subsets including human pan-T cells, CD4+ T cells, CD8+ T cells, B cells, and NK cells, with 49 recombinant chemokines using a singular technique, and standardized conditions resulting in a dataset representing 238 assays. We then outline a groundwork computational model that can simulate cellular migration in the tumor microenvironment in response to a chemoattractant gradient created from stromal, lymphoid, or antigen presenting cell interactions. The resulting model can then be parameterized with standardized data, such as the dataset presented here, and demonstrates how a computational approach can help elucidate developing ELNs and their impact on tumor progression.

Isolation of Lymphocyte Populations. Spleens and the superficial inguinal, axillary, lateral axillary, mesenteric and cervical lymph nodes of mice were removed under sterile conditions and mechanically dissociated to prepare single-cell suspension with a 100-µm nylon mesh. Erythrocytes were lysed with RBC lysing buffer (0.15 M NH 4 Cl, 1 mM KHCO 3 , and 0.1 mM EDTA in sterile water). After incubation for 1 minute, cells were washed in 1x PBS and used for further experiments.
Pan T, CD4 + T, CD8 + T, B and NK cells were then isolated from splenocyte and lymph node cell suspensions using negative immunomagnetic isolation kits as per the manufacturer's instructions: Pan T Cell isolation kit II, CD4 + T Cell isolation kit II, CD8α + T Cell isolation kit II, B Cell isolation kit and NK isolation kit for mouse (Miltenyi Biotech, Auburn, CA). Isolations were performed using a magnetic cell sorter (autoMACS) (Miltenyi Biotech, Auburn, CA). Pan T and B cell isolation enrichments were greater than 95%, CD4 + T and CD8 + T cells were greater than 90% and NK cells were greater than 80% as determined by flow cytometric analysis. See supplemental Fig. 1.
Microchemotaxis Assay. Each chemokine was added to the lower chamber of a 24-well transwell (Costar, Cambridge, MA) at indicated concentrations in 600 μL of CM, and incubated at 37 °C in a humidified incubator with 5% CO 2 for 30 minutes. Respective lymphocyte populations were resuspended at 1 × 10 7 cells/mL in CM, and 100 μL (1.0 × 10 6 cells) were seeded into 6.5-mm 24-well transwell inserts (Costar, Cambridge, MA), and allowed to incubate at 37 °C for 10 minutes. After pre-incubations, the upper chambers were moved into lower wells to start assay. To determine the actual input amount, 100 μL of cell suspension was added to 500 μL of CM and put into the lower wells directly and incubated without upper chambers. After 3 hour incubation at 37 °C, the upper chambers were removed and the migrating cells were retrieved. Cells were stained with the appropriate antibodies and added 2 × 10 5 polystyrene beads (Bangs Laboratories, Fishers, IN). Twenty thousand beads were then counted flow cytometrically. Baseline chemokinesis was determined by the fraction of cells moving into the lower chamber in control wells containing no chemotactic agent, and is expressed as a percent of seeded cells. The chemotactic index (CI) is calculated as the fraction of cells migrating in response to the chemotactic agent normalized to the chemokinesis background such that CI = (fraction of cells migrating to condition)/(fraction of cells migrating to control media). Transwell migration through 3 μm, 5 μm, and 8 μm pore sizes was compared (Supplemental Fig. 2), and the 5 μm pore size was selected for use. Data reported are representative of at least two experiments.
Enzyme-linked Immunosorbent Assay (ELISA). ELISA was carried out to evaluate T and B cell activation (Supplemental Fig. 3D). After supernatants were harvested at day 5 in the T cell activation assay or at day 3 in the B cell activation assay, IFNγ and IgM were measured, respectively. Mouse IFNγ ELISA Kit II (BD Biosciences) Figure 1. Resting NK cells display higher chemokinesis than other lymphocyte populations. Chemokinesis is reported as the baseline random migration into the lower chamber of a control transwell assay containing no chemoattractants. ****p < 0.0001. and Mouse IgG ELISA Quantitation Set (Bethyl laboratories, Montgomery, TX) were used in this experiment according to the manufacturer's instructions.
Statistical Analysis. Unpaired student t-test was performed for each chemokine concentration versus control wells containing no chemokine using GraphPad Prism 5.04 software (GraphPad Software, La Jolla, CA). Microchemotaxis assays were performed on individually isolated/activated wells for each cell population (n = 4).

Mathematical Model.
A cellular automaton model was developed to investigate the effect of ELNs on the organization and effectiveness of the immune response to a nearby tumor. The mathematical model is simulated in two dimensions and contains a tumor, a patch of fixed RFC, and four types of motile cells that move in a continuous domain (i.e., off-lattice): inactive T cells, activated T cells, resting antigen presenting cells (APC OFF ), and activated marginal zone antigen presenting cells (APC M ). The tumor is represented by a circle of radius R, where  T cell response to CCL19, CCL21, and CXCL10 at the indicated concentration ranges. (C) Resting enriched CD8 + T cell response to CCL19, CCL21, CXCL10, and CXCL11 at the indicated concentration ranges. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Scientific REPORTS | 7: 15996 | DOI:10.1038/s41598-017-15924-2 we assume that the radius grows linearly in time in the absence of any immune response, and that the T-cell response causes a reduction in growth rate proportional to the number of tumor-infiltrating lymphocytes (L). The equation for tumor growth per time step of the simulation is given by: where g is the innate tumor growth rate, k is the killing rate of TILs, and Δt is the time step used in the simulation. Assuming a disease state in which anti-tumor cytotoxicity can occur, increasing numbers of TILs would slow the growth and potentially cause the tumor to regress if kL > g. APC OFF are introduced into the tumor area at a constant rate. The new cells are randomly positioned in an annulus centered around the tumor, with inner radius of (R − 500 µm) and outer radius of (R + 1500 µm). Once in the environment, they move with a combination of an unbiased random walk and a directed walk towards the tumor center, following tumor-produced inflammatory chemokine gradients. When an APC OFF is inside the tumor, the antigen collection process is simulated by converting the cell to an APC M , which can present antigen to the T cells. These APC M are also motile, using an unbiased random walk in combination with a directed walk away from the tumor center, to represent the seeking of T cells and/or vasculature for subsequent antigen presentation.
Inactive T cells are added to the simulation with a constant rate, using the same annulus as for APC OFF cells. These cells have an unbiased random walk. When an inactive T cell encounters an APC M (based on the two cells being in proximity to each other at a given time step), the APC M activates the T cell. The APC M can activate a number of T cells before it is removed from the simulation. Active T cells follow a biased random walk in the direction of the tumor. When they cross the tumor boundary as determined by R, they become TIL, which affects the tumor killing rate as shown in Eq. 1. They remain TIL for a set period of time, after which they are removed from the simulation, as expired T cells have no further effect on the results.
To model the effect of ELNs, we randomly place a fixed number of RFC in a circular patch situated away from the tumor. Each RFC begins in the "off state" with no secretion of chemokines, and therefore no influence on the movement of any of the APC or T cells. However, when an APC M migrates and comes into "contact" with an RFC (based on the two cells' proximity to each other), it will activate the RFC, which will in turn start to produce a chemokine gradient. For simplicity here, we model the chemokine gradient produced by an activated RFC as a two-dimensional Gaussian function centered at the cell. The total chemokine gradient for multiple activated RFCs is the sum of these Gaussians. As more RFC are activated, the chemokine signal becomes stronger. This gradient affects two of the motile cells, APC M and inactive T cells. For both of these cell types, when they detect the presence of chemokine gradient produced by RFC, movement is biased toward the ELNs patch. The bias scales with the strength of the gradient. We vary the number of RFC in different simulations to investigate this effect.
The random walks in the model were implemented as follows: the motile cells in the simulation are off-lattice and therefore can take any value for their position. The maximum travel distance for a cell in a time step is its kinetic speed multiplied by the time step. The algorithm selects a random distance between 0 and this maximal distance, and moves the cell in the direction of a randomly chosen angle. When there is an additional chemotactic term, a second random distance is chosen based on the maximum tactic speed, but the angle is fixed, directed towards the chemokine source (ascending the gradient). In this model, we do not consider spatial occupancy, so cells can get arbitrarily close to each other. To summarize the model, there are two chemokine gradients, one from the tumor and one from the RFC M . APC OFF and active T cells move towards the tumor, and APC M and inactive T cells move away from the tumor. When there are activated RFC M , APCM and inactive T cells are both drawn to the ELNs. The simulation is initialized with a small tumor, a patch where RFC OFF cells are placed, and 30 APC OFF and inactive T cells. Parameters used in the model are shown in Table 1, with those labeled as 'model specific' refer to parameters that are highly variable in practice (e.g., tumor size, ELN-tumor distance, influx rate of immune cells); values were chosen that were in line with biologically reasonable ranges.

Results
Resting pan T cells, B cells, CD4 + T cells, CD8 + T cells, or NK cells were immunomagnetically isolated and chemotaxis assays were performed for the 48 chemokines listed in Table 2 at three different concentrations: 10 ng/ mL, 100 ng/mL, and 1,000 ng/mL. Non-directional random motility, or chemokinesis, was assessed for each cell population using control wells containing no chemokine gradient. Resting pan T cells, CD4 + T cells, CD8 + T cells, and B cells all displayed relativity low levels of chemokinesis while NK cells displayed significantly higher non-directional movement (Fig. 1). Directional movement in response to a chemokine concentration gradient was normalized to baseline chemokinesis for each respective cell type, and is reported as a chemotactic index (CI) for each condition. CI was calculated as the fraction of cells that migrated into the lower chamber from the total number of cells seeded in the top chamber for each condition, and then divided by the fraction of cells that migrated into the lower chamber randomly in control wells (i.e. by chemokinesis).
When activated, lymphocytes become more motile and often down-regulate expression of chemokine receptors recognizing secondary lymphoid tissue-homing chemokines such as CCR7. This aids in emigration from the follicle and secondary lymphoid tissues 34 . Here, activation of pan T cells, enriched CD4 + T cells, and enriched CD8 + T cells with plate-bound anti-CD3 resulted in significantly increased chemokinesis. Similarly, but to a lesser extent, activation of B cells with LPS significantly increased chemokinesis (Fig. 5A). Despite higher baseline motility, activated T cells loose chemoattraction to lymphoid-homing chemokines as demonstrated by reduced CI toward CCL19 and CCL21 for pan T cell populations (Fig. 5B,C), and reduced chemoattraction of enriched CD4 + and CD8 + T cells toward CCL21 (Fig. 5D,E). In contrast, B cells displayed significantly increased chemoattraction toward CCL21 and CXCL13 upon activation with LPS (Fig. 5F,G). This likely reflects the divergent roles of B cells and T cells following activation in secondary lymphoid tissues with T cells often emigrating T-cell zones and recirculating into the periphery and with B cells clustering in follicles.
The standardized set of chemotaxis assays presented here demonstrate that chemokines differentially affect random and directed motion of lymphocytes. To investigate the implications of this regarding immune response and subsequent tumor growth depending on the presence of the ELNs we developed a phenomenological mathematical model of the tumor and surrounding environment that is able to make use of chemoattraction data an simulate different scenarios. Figure 6 shows snapshots of the model at 14 days of simulation, under three different conditions for the RFC: (a) no RFC, (b) 5 RFC, and (c) 30 RFC. Time course plots for these and several other simulations with different RFC counts are shown in Fig. 7, in the left panel. The results suggest that the presence of the RFC has an impact on the immune response to the tumor, both in terms of T-cell activation and tumor reduction. In all cases, the tumor grows for an initial period since the immune system starts in an inactive state. Eventually, APCs present tumor antigens to inactive T cells, causing activation, migration, and accumulation of TILs inside the tumor and eventual tumor regression. The presence of higher numbers of RFC causes a faster immune response to occur: without an ELNs, the peak tumor size occurs at day 9; while for simulations with 100 RFC the peak tumor size is around 4 days, suggesting that the migratory organization of lymphocytes provided by the ELNs significantly accelerates the anti-tumor immune response. In the case with 100 RFC, the tumor also reaches a maximum size that is about 25% smaller than the simulation with no RFC. This advantage persists as tumor growth continues to decline. The right panel of Fig. 7 compares tumor sizes after 30 days of simulation. Of interest, having a small number of RFC cells (5-10) does not produce an immune response as effective as the case with no RFC or a high number of RFC. This nonlinear result is because the weak chemokine signal of the ELNs draws some nearby cells away from the tumor, but a lack of immune cell density limits the effectiveness of rapid activation. In essence, there is a need for critical mass in the ELNs to achieve efficient cell attraction and activation. This can be seen in Fig. 6B, where the dendritic cells are accumulated in the ELNs area, but inactive T cells remain predominantly in the tumor microenvironment due to the weak signaling from the ELNs patch.

Discussion
We had earlier identified a unique tumor-derived, 12-chemokine gene expression signature that could accurately predict the degree and type of lymphoid infiltrate, organized remarkably as ELNs that comprise -by immunohistochemistry staining -prominent B cell follicles, T cell marginal zones comprised of both CD4 + and CD8 + cell subsets, and associated follicular dendritic cells 23 . Of importance, there was a highly significant and consistent association between a marked increase in overall patient survival, the value of the mean score of this gene expression signature, and the presence of ELNs in stage IV (non-locoregional) melanoma, colorectal cancer, and stage IV bladder cancer. However, we found that the majority of human solid tumors lacked the presence of these particular ELNs, and patients harboring these 'poorly-immunogenic' tumors have had uniformly poor prognosis (i.e. reduced overall survival). Thus, there is a clear unmet medical need to 'reverse' this tumor microenvironment (that lacks these particular ELNs) by manufacturing de novo 'designer' ELNs.
Lymphocyte recruitment to the tumor microenvironment represents an attractive target to enhance anti-tumor immunity. As we continue to study lymphocyte trafficking, it is increasingly necessary to view chemokine-mediated trafficking as networks of chemokines and chemokine receptors working in concert as opposed to individual chemoattractive axes. Understanding how chemokine receptor expression and lymphocyte responsiveness changes during the course of maturation, activation, and effector response will greatly inform the development of preclinical animal and mathematical models. The chemokine database described herein may serve as a normalized resource to parameterize such models. Given a direct link between gene expression of certain chemokines within human solid tumor microenvironments, the presence of tumor-localized ELNs, and addressing a clinical unmet need, we embarked on developing a standardized database of chemoattractive potentials of immune cell subsets for a broad range of chemokines to identify candidates for future engineering of ELNs. To approach this goal, we developed a strategy that combined the use of defined, recombinant chemokines, highly enriched resting and activated immune cell subsets, with standardized microchemotaxis assays to provide potential leads to further test in mathematical models. In this regard, we are developing both mathematical and pre-clinical animal models of ELNs formation in which multiple elements are being interrogated. The inclusion of lymph node-derived primary cellular components, which normally provide chemotactic and homeostatic queues in conventional lymph nodes, are being genetically modified to express selected chemotactic (e.g., see Table 3) and lymphoid neogenesis-related genes to enhance ELN formation. These modified cell lines are being combined with tumor antigen-pulsed dendritic cells and then incorporated in biocompatible scaffold materials and administered to tumor-bearing mice as injectable or implantable matrices 35 . These matrices may provide two potential benefits when delivered to tumor-bearing hosts. First, they may serve as model systems to better understand the factors governing the formation and/or maintenance of ELNs with anti-tumor reactivity. Second, these matrix-based systems may function as a therapeutic platform by delivering, stimulating and expanding transplanted lymphocytes and/or dendritic cells. In addition, the inert nature of bio-scaffolds also allows for the implementation of microparticle or nanoparticle constructs for controlled release of soluble factors. Such measures can provide sustained environmental queues to augment antigen-presenting cell or lymphocyte longevity, maturation, and activation. The initial mathematical modeling presented herein will be further developed to include these parameters.
In regard to translating in vitro migration assays to in silico models, care must be taken to differentiate between non-directional movement known as chemokinesis and directional chemotaxis toward a chemokine gradient. CI, as it is calculated here, normalizes transwell migration relative to baseline chemokinesis. This denominator is greatly increased after lymphocyte activation. Without normalizing to baseline movement, the percentage of cells migrating may appear very high when the majority of migration can be accounted for by  random motility. Because this database was constructed using a single methodology, time point, and concentration range, CI can be compared across chemokines. For example, CCL19 and CCL21 are both well studied T cell chemoattractants; however, here CCL21 appears to be a 50% stronger chemoattractant compared to an equimolar amount of CCL19. In addition, this database serves to point out suboptimal chemoattractants whose combined cumulative effects may play an important role for in silico models such as the broad responsiveness observed for NK cells.
As a first step toward integrating this type of chemotactic data into an in silico model of ELNs formation, we developed a simple phenomenological mathematical model that predicts the chemokine gradient created as a result of lymphocyte, APC, and stromal interactions in the tumor microenvironment. Here, we consider secretion of a general chemoattractant for responding lymphocytes, and can serve as the basic groundwork for simulating multiple chemokine/chemoattraction axes involving more detailed cellular phenotypes. Further investigation of relevant chemokine relations between APC, RFC, and T cells would lead to a more mechanistic implementation of the model, which could then inform the design of future in vivo studies. The generalized gradient presented here can be replaced by specific chemokines observed in the microenvironment of tumor samples, or with chemokines that are known to be secreted by activated DC and/or stroma. In this way, multiple gradients can be modeled in concert. In addition, each chemokine gradient can be further developed beyond a two-dimensional Gaussian distribution taking into account extracellular matrix components that can bind to and slow the degradation of chemokines, better representing lymphocyte conduit systems observed in follicular structures in vivo 28 . So too can the responsiveness of activated cell populations be adjusted based on lymphocyte activation, and on negative feedback mediated by shingosine-1-phosphate receptor (S1P 1 ) activity [36][37][38] . Importantly, this model outlines the framework to simulate chemotactic movement of infiltrating populations in response to chemokine gradients. This alone inadequately represents the physiologic tumor microenvironment, and in particular does not address immune cell polarity or immunosuppression. In addition to representing more defined subpopulations of immune cells, other soluble factors governing tumor-induced immune suppression will likely need to be added to fully capture the dynamic interplay of pro and anti-tumor elements. In particular, TFG-β, IL-10, and regulatory T cell populations are prevalent components of the tumor microenvironment known to suppress anti-tumor immunity 39 . Concurrently, the production of type I-polarizing factors such as IFN-γ and IL-12 can also be added as a counterbalance to these type II-polarizing factors 40 . Furthermore, it is now appreciated that the tumor as a system is heterogeneous, composed of numerous dominant and subdominant clones with related but different genetic lesions and mutational loads 46 . This diversity results in compartmentalized spatial fields across the tumor bead, each with its own variation of environmental queues, grow rates, and immune involvement 47 . It is reasonable to expect that variations in the type and amount of soluble factors released can also vary across the tumor bed in accordance with clonal diversity. Such variations would directly impact the localization and function of infiltrating immune cells, and may be required to fully model the complex tumor physiology.
The addition of these factors in future models would result in a system whose outcome is not guaranteed, but would rather depend on the balance of pro and anti-tumor forces recruited or induced in the tumor microenvironment. The complexity such models will increase exponentially with the addition of individual chemokine gradients, further subdivision of inflammatory or suppressive cell types, the addition of polarizing soluble factors, and representation of clonal heterogeneity, necessitating the use of standardized datasets such as is presented here. The model introduced here, though simplified, clearly demonstrates the potential benefit of ELNs and suggests that stromal/APC interaction is paramount for effector lymphocyte organization and response. Continued development of integrated ELNs formation models will greatly inform potential points of therapeutic intervention, and generate novel hypotheses regarding anti-tumor immunity.  Table 3. Chemotactic Index (CI) Summary for the 12-chemokine GES on Resting Lymphocytes. "−" No significant difference in CI at any chemokine dose. "−/+" p < 0.05 for at least one chemokine dose, but CI less than 5 for at least one chemokine dose. "+" p < 0.05 for at least one chemokine dose, and CI over 5 for at least one chemokine dose lg. ND Not determined.