Autologous gradient formation under differential interstitial fluid flow environments

: Fluid flow and chemokine gradients play a large part in not only regulating homeostatic processes in the brain, but also in pathologic conditions by directing cell migration. Tumor cells in particular are superior at invading into the brain resulting in tumor recurrence. One mechanism that governs cellular invasion is autologous chemotaxis, whereby pericellular chemokine gradients form due to interstitial fluid flow (IFF) leading cells to migrate up the gradient. Glioma cells have been shown to specifically use CXCL12 to increase their invasion under heightened interstitial flow. Computational modeling of this gradient offers better insight into the extent of its development around single cells, yet very few conditions have been modelled. In this paper, a computational model is developed to investigate how a CXCL12 gradient may form around a tumor cell and what conditions are necessary to affect its formation. Through finite element analysis using COMSOL and coupled convection-diffusion/mass transport equations, we show that velocity (IFF magnitude) has the largest parametric effect on gradient formation, multidirectional fluid flow causes gradient formation in the direction of the resultant which is governed by IFF magnitude, common treatments and flow patterns have a spatiotemporal effect on pericellular gradients, exogenous background concentrations can abrogate the autologous effect depending on how close the cell is to the source, that there is a minimal distance away from the tumor border required for a single cell to establish an autologous gradient, and finally that the development of a gradient formation is highly dependent on specific cell morphology.


Introduction
Within the brain many dynamic processes take place on a regular basis in order to maintain homeostasis. These include electrical gradients being induced and sent out to the rest of the body for communication and motor function, fluid exchange in order to keep the correct amount of ions, nutrients, and oxygen delivered to the tissue, and biophysical forces acting on the tissue to direct organization and repair [1]. On a smaller scale, these processes are driven by interactions at the cellular level and governed by signaling cascades, cell-to-cell interactions, cell-to-environment interactions, chemical gradients, and individual cellular phenotypes/functions. At this level, biophysical forces such as fluid flow have a direct impact what a cell feels and how it reacts, for example flow-mediated shear stress on the outside of a cell leading to cytoskeletal remodeling and cell movement [2], [3]. Another example is flow acting as a convective force to create gradients of chemokines which can help pattern organs such as during neonatal development [4], [5], recruit immune cells to distant sites of inflammation [6], or even play a role in tumor invasion and development [7], [8]. To complicate matters, this flow is driven by multiple means such as respiration and blood pressure, through multiple structural pathways such as blood vessels, lymphatics, and ventricles, and is constantly changing in magnitude [9]. In the brain, the fluid that cells in the parenchyma encounter is interstitial fluid. This is similar in property to cerebrospinal fluid and is a combination of fluid coming from blood vessels and cerebrospinal fluid from the ventricles into the Virchow-Robin space [10]. This fluid travels through the extracellular matrix (ECM) and around cells, driven by the pressure differential from the vasculature. More recently, the glymphatic pathway has been characterized as a conduit by which the interstitial fluid flows in/around arteries, through the extracellular space, and ultimately is drained by the cervical lymphatics as it travels around veins to leave the cortex [11]. In this context, we can think of interstitial fluid as having a defined, or at least somewhat organized, direction of flow governed by anatomical structures that also directly interacts with cells.
Interstitial fluid flow is especially of interest because it has been shown to change drastically with pathological states such as cancer, hydrocephalus, and Alzheimer's disease [10], [12]. In particular, our lab and others have identified interstitial fluid flow as a mediator of cancer cell invasion in brain and breast cancer cell lines [7], [13], [14]. We have also mapped velocities based on human MRI around gliomas and shown that magnitudes and flow direction are heterogeneous with average speeds of 3um/s [9]. During tumor development, intratumoral pressure increases while pressure in the surrounding tissue stays at normal levels [15]. Specifically, the pressure stays mostly constant in the tumor interior and drops quickly at the tumor border. This pressure is a result of leaky vasculature caused by breakdown of the blood brain barrier and remodeling of the ECM by the tumor. The resulting pressure differential causes interstitial fluid flow from the tumor into the surrounding space, with implications for tumor growth and development [7], [16], [17]. This is especially important in glioma as recurrence is caused by cells that have invaded past the tumor border and evaded removal by surgical resection because they are outside of the surgical margin. This brings up important questions as to what is causing this migration, especially in regards to fluid flow effects. There are currently two main hypotheses by which flow as a biophysical force is mediating this invasion: mechanotransduction and autologous chemotaxis. The model in this paper is focused on autologous chemotaxis and how chemokines form under different flow conditions.
Chemokines are chemoattractant cytokines found ubiquitously throughout the body. They are involved with processes such as cell proliferation, cell migration, and the inflammatory response. In the context of cancer, chemokines have been found to promote invasion and tumorigenesis [18]. There are over 40 main chemokines not including their receptors [19] that cover a range of functions and are found in many different types of cancer. Some of the more well-documented include CCL21 [20] in breast cancer, CXCL12/CXCR4 in brain cancer [19], and CCL2 in prostate cancer [21]. These chemokines are also not exclusive to any one type of cancer. For example, CXCL12 induces EMT in colorectal and pancreatic cancer cells and has also been implicated in neoangiogenesis or angiogenic activity through PI3K pathway in glioblastoma [19]. Because of its high levels in glioma and the many roles that it plays with promoting glioma progression, we chose to focus on CXCL12 and its receptor CXCR4. It is also important to realize the way in which this axis interacts with the cell and surrounding matrix. As CXCL12 is released into the environment, it not only binds to specific sites on the cell surface but also molecules within the extracellular matrix. One of the main molecules is a glycosaminoglycan (GAG) called heparan sulfate (HS) that allows the CXCL12 to bind directly and plays an important role in regulating chemokine gradients within the microenvironment [22]. In our model we are specifically interested in the bound CXCL12 gradients that form as a result of the relationship with HS as this ultimately dictates cell movement that occurs from the CXCL12/CXCR4 axis. In addition to secreting chemokines, glioma cells also secrete proteases at higher levels than normal cells [23], [24]. These substances, such as MMP9, can break down the extracellular matrix, both creating spaces for the cells to migrate as well as releasing bound molecules such as the CXCL12.
In 2006, Fleury et al. published a computational paper on the theory of autologous chemotaxis as a result of fluid flow [25]. They investigated autologous chemotaxis in the context of CCL21 and CCR7 through a 2D finite element model. They incorporated unidirectional fluid flow across a cell (represented by a circle) with secretion of a chemokine and a protease from the cell, while looking at the resulting pericellular gradient that formed under different reaction conditions and velocities. Their model found that gradient formation is dependent on interstitial fluid flow velocity, increasing drastically from just 0 um/s to 6 um/s flow magnitude. This model was expanded by Evje [26] in 2018 to include multiphase components (treating the fluid and the cancer cells as two distinct phases with their own equations) and adding a stress term to describe the force the fluid exerts on the cell to try to "push" it downstream.
In 2007, Shields et al. [20] published a seminal paper on autologous chemotaxis, with the hypothesis that cells migrate along autocrine chemokine gradients formed by fluid flow from the lymphatics in breast cancer/melanoma. They specifically looked at CCL21 and its receptor CCR7 and found that chemotaxis toward lymphatic endothelial cells was mediated by CCR7 and that flow induced a higher % migration dependent on CCR7 but independent of LECs, suggesting that a mechanism such as autologous chemotaxis is present.
There have been other papers published recently that investigate gradient formation of certain chemokines in the context of inflammation, the immune system, and other diseases [27], [28] but none that investigate the role of CXCL12 in autologous chemotaxis. We build on the foundation of previous models here by bringing it into the 3D realm and considering effects such as cell shape and time on gradient formation, which has not been done before to our knowledge. We also extend the application to glioma and the CXCL12/CXCR4 chemokine axis as a way to describe how tumor cells might be migrating in the brain due to changes in interstitial fluid from tumor growth and drainage via the glymphatic system.

Materials and Methods
Equations and implementation of the model Our model focuses on a control volume around a single cell and looks at gradient formation under different parameter changes, flow directions and magnitudes, cell morphology, transient functions, background concentrations, and distances from the tumor border. COMSOL is used as the software package in which to set up the finite element analysis and run the model. The model is based on modified Navier Stokes (convectiondiffusion equation): where u is the velocity, v is the dynamic viscosity, p is the pressure, and g is gravity. The first term describes the change in velocity over time, the second term the convection component, and the third term the diffusion component. The terms on the right hand side make up the internal and external forces on the system, with a term for the pressure gradient and a term for gravity. This can be simplified by the assumption of stationary (no time dependence), creeping (convection small compared to viscous or diffusion), and incompressible flow (ρ is constant) to the Stokes equation: from this the Brinkman form of Darcy's law can be derived by assuming viscous force is linear with velocity and adding the Brinkman term to account for transitional flow. This equation describes fluid flow through a porous media (such as interstitial fluid flow through the ECM): 4 where beta is a correction term to account for flow through porous media, q is the flux (flow per area), k is permeability, µ is the viscosity, u is the velocity, and p is pressure. This equation is coupled with the mass transport equation which describes transport of a chemical species through convection and diffusion along with any reactions taking place: where c is the concentration and D is the diffusion coefficient. The first term describes the change in concentration over time, the second term the convective component of concentration change, the third term the diffusive component of concentration change, and the last term the chemical reactions that take place to either generate or degrade the chemical species. The reaction equations used in this model describe the binding kinetics of CXCL12 to the matrix. As the CXCL12 is secreted from the cell, it preferentially binds to sites on the extracellular matrix, such as heparan sulfate [22]. This bound CXCL12 is what will ultimately interact with the CXCR4 receptor on the cell to cause migration or other cell activity and is the main concentration that we are concerned with in the model. The CXCL12 binding is governed by the kon and koff values, experimentally determined by Munson et al [17]. In addition, the cell secretes a protease, such as MMP9, that acts to cleave the bound CXCL12 and release it back into the environment as free CXCL12. The protease is an enzyme in the model, meaning it does not get altered by its proteolysis function and so its reaction is zero. The equations are shown below: Protease: = 0 where [cxcl12] is the concentration of the free chemokine, [hs_cxcl12] is the concentration of the bound chemokine, [HS] is the heparan sulfate concentration, and [protease] is the concentration of the protease. The solution for the Brinkman equation (3) is first generated and the corresponding velocity term is used by the mass transport equation (4) to solve for the concentration at each element in the model. Mass transfer coefficients for the chemokine and protease were calculated via correlation 'forced convection around a sphere' [29]: where km is the mass transfer coefficient, D is the diffusion coefficient, Dsphere is the diameter of the sphere, ρ is the fluid density, µ is the fluid viscosity, and V is the average fluid velocity.
In order to quantify model results, % concentration is calculated as: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 November 2021 where Ci is the concentration taken at the respective point on the surface of the sphere, either upstream or downstream of flow. All parameters and their values for the model can be found in Table 1.

Conditions and Quantitative
Values used in the model Two COMSOL modules are used in this model: 'fluid flow' and 'chemical species transport'. The 'Brinkman Equations' physics and 'Transport of Diluted Species in Porous Media' physics are then used to set up the specific equations for the model. Two domains are present in this system -one for the cell, modeled as a sphere, and one for the control volume, modeled as a rectangular prism. COMSOL Multiphysics is used to establish the relationship between the brinkman equation and the mass transport of the CXCL12. The HS is modeled in a separate physics for transport where there is no convection or diffusion (the HS is present through the ECM as a structure that does no move). There are many boundary conditions prescribed throughout the model, with each figure having a slightly different configuration. For the baseline case, a no slip boundary condition is given for the fluid flow around the outside of the cell, thus there are viscous effects at this interface (u = 0 m/s at this boundary). An open boundary condition is used at the walls of the control volume for the fluid flow equations as this allows movement of fluid into and out of the domain. An inlet and outlet Dirichlet boundary condition are specified at either end of the control volume (Supplemental Figure S1) with the inlet being a uniform velocity (u = 1E-5 m/s) and the outlet being a uniform pressure (p = 0 Pa). From the mass transport side, a Neumann boundary condition is applied around the outside of the cell representing secretion of the chemokine CXCL12 and the protease MMP9. An open boundary condition is specified for all of the outer walls of the control volume where convective inflow and outflow of the chemical species can occur. This represents an open space where the chemokine can disperse as if the surrounding space is a large volume, which is applicable for our system. Reactions (Equations 5, 6 & 7) are prescribed to take place within the control volume. Initial values for the chemokine and protease were set to 100 nM and 1 nM, respectively. Geometry of the figure is set up to include a rectangular prism 100 µm x 50 µm x 50 µm (LxWxH) corresponding to a control volume encompassing a cell, which is modeled as a sphere with diameter of 10 µm. The sphere is placed 25 µm from the inlet boundary. Meshing was done through the COMSOL software to apply a free tetrahedral mesh to the geometry. A predefined mesh calibrated for fluid dynamics was used with a 'finer' element size applied for each domain based on a mesh refinement analysis (Supplemental figure S2). For specific model conditions applicable to individual figures, refer to the supplementary materials.

Assumptions of the model
The model presented here has inherent assumptions that need to be addressed. Interstitial fluid is treated as a Newtonian fluid in this model and the flow through the system is assumed to be laminar, viscosity-dominated, and incompressible (Stokes flow). We also model the cell as continuously secreting both the chemokine and the protease. Under normal circumstances this secretion may fluctuate based on factors in the environment such as shear stress exerted on the cell or individual cell phenotypes. The reaction equations represent simplified relationships between the CXCL12, protease, and HS. While values from the literature were used where possible, some terms such as the concentration of the HS and the krel term were guessed. These parameters have not been well-studied and should be experimentally determined for exact values in the future, especially because a sensitivity analysis reveals that krel accounts for large fluctuations in model outcome (Supplemental Figure S3). In addition, we only model one protease, MMP9, and one chemokine, CXCL12. There are many different matrix metalloproteinases that may be at work, however, to impact the gradient formation based on ECM degradation. Lastly, the control volume that we have chosen was assigned uniform parameters for porosity and permeability which do not change in relation to the protease in contrast to what would be seen with ECM remodeling in vitro or in vivo. We do not include a term for the protease degradation and we do not take into account any interaction between the cell and the chemokine such as recycling effects or receptor binding. More detailed information specific to each manipulation can be found in the Supplemental Methods.

Effects of transport parameter changes on pericellular gradients
To begin, we probed the effects of changing biotransport parameters on pericellular gradients. From parametric sweeps of each individual parameter, it is clear that three in particular stand out as affecting gradient formation the most (Supplemental Figure S4): diffusion coefficient ( Figure 1A), reaction coefficients for CXCL12 binding to HS to form bound CXCL12 (1C), and velocity of flow (1E). The values used were from a physiological range pulled from the literature ( Table 1). The resulting concentration gradients at steady state were plotted as a function of the distance from the cell. We found that for diffusion coefficient (1B) the bound CXCL12 profile changed more in magnitude than %concentration showing that while diffusion coefficient might affect the overall concentration of the free and bound CXCL12, it does not significantly impact the gradient that forms around the cell. A smaller diffusion coefficient of CXCL12 results in a greater accumulation of the CXCL12 and therefore more bound CXCL12 and a more distinct pericellular gradient. Interestingly, a smaller diffusion coefficient of protease results in a lower bound CXCL12 concentration and less downstream gradient formation. This makes sense if we consider that a more concentrated protease (i.e. lower diffusion) will cleave more of the bound CXCL12, disrupting the gradient formation. For reaction coefficient (1D) the binding and cleaving rates of the CXCL12 to the HS impacts the resulting pericellular gradient immensely. If the release rate is increased we see a lesser concentration of bound CXCL12 and eventually a disruption of the gradient as the bound CXCL12 is mostly cleaved to become free CXCL12. Changing the kon or koff values impacts the magnitude of the bound CXCL12 -increasing kon will increase the bound CXCL12 and increasing koff will decrease the bound CXCL12 as expected. Lastly, for velocity we saw that higher velocities led to steeper differences of concentration upstream and downstream of the cell (1F). As velocity is increased, so too is the %concentration as the chemokine is affected by the convection of fluid flow. The difference across the cell can be used to calculate a concentration gradient (1G), which can be used to simply describe the steepness of the gradient across the cell and thus the gradient that a cell may feel. This value is denoted as %concentration and will be used throughout the manuscript. As expected, we see a significant correlation (r=0.991, p<0.05) between the velocity and the %concentration (1H).
Velocity has the largest impact on % concentration of bound CXCL12 with 0 m/s corresponding to about 0% all the way to 100 µm/s causing a gradient of 46%. Interestingly, it's not until 10 µm/s that we see a large increase in gradient formation occurring (>5%). This is apparent in Figure 1E where we can see the concentration gradient profiles and the increase in the formation of bound CXCL12 skewed to one side of the cell with increasing velocity. At 0.1 and 1 µm/s the % concentration is 1% or lower. The results found here are in agreement with those found by Fleury et al.

Directionality of flow alters gradient
Magnitude of interstitial fluid velocity corresponded with increasing concentration gradients across the cell. However, in vivo, we see not only a range of interstitial flow velocity magnitudes, but also direction of the vectors. For example, in mice, we have used dynamic contrast enhanced MRI and mass transport analysis to map interstitial fluid flow velocities across tumors and interstitial space in the brain (Figure 2A). As the tumor progresses, fluid flow becomes heterogeneous especially around the tumor border where invading cells are subjected to its effects. These flow patterns calculated from MRI can be overlaid on histological sections to examine cellular invasion within the brain tissue beyond the tumor border. The resolution of the velocity vector field is 104.2 μm/pixel. We identified and modeled two individual flow patterns (2A.1 and 2A.2): one where the flow vectors came from opposite sides of the cell and one where the flow was split into an xand y-component with higher flow in the y-direction (2B, C). Introducing flow direction from either side of the cell lead to an even distribution of bound CXCL12, with no chemokine gradient formation. Interestingly, while no gradient developed, we did see an induced surface pressure across the cell which might have implications in other cell migration mechanics. When fluid is induced in multiple directions, the formation of the gradient follows the resultant of the flow direction (2C). We have observed areas of higher invasion in regions of outward fluid flow from the tumor bulk [17], so the development of the gradient in the direction of flow seems to support the idea of autologous chemotaxis. We have also modeled different flow arrangements and have seen this same effect of multidirectional flow skewing the gradient toward the resultant flow direction (Supplemental Figure S5).

Temporal fluctuations in velocity yield variable gradients
After exploring some of the effects of IFF velocity magnitude and direction, we wanted to understand the time-sensitive effects on gradient formation. Transient solutions were first solved at different time scales in order to investigate the rate at which the gradient formation could actually be observed. Time scales at minutes, seconds, and centiseconds were computed. As shown in Supplemental Figure S6, the gradient can be observed to develop on the centisecond time scale. This seems to be within ranges that have been observed for ligand binding [30]. We next wanted to test the effect that varying flow magnitudes has on the developing cell gradient. We settled on three different conditions that would be physiologically relevant to brain cancer: convection enhanced delivery (CED), surgical resection, and pulsatile flow from cerebrospinal fluid (CSF). CED is a treatment method that allows for the bypassing of the blood brain barrier to deliver therapeutics via catheter directly to the area of the tumor. This method introduces a volume of fluid, which actually increases the interstitial fluid flow throughout the parenchyma. This was modeled as a ramp function with increasing IFF velocity over a short period ( Figure 3A) which corresponds to an increasing %concentration that looks almost logarithmic in nature. This is unsurprising based on the results from Figure 1, but we can see that at some point, the gradient starts to drop off as it approaches a potential maximum threshold around 80%. Surgical resection causes a large drop in pressure based on experimental results [31] which we modeled here as a step function in the IFF velocity (3B). The initial step up shows IFF increase due to tumorigenesis with a steep drop off when resection happens. We see the corresponding rapid increase in bound CXCL12 gradient as the velocity increases, but then almost half a reduction in %concentration as the velocity returns to a more normal value (which would be associated with the resection). Interestingly, the %concentration does not return to normal but instead retains an increased value from baseline. During normal function, the brain has oscillatory flow due to CSF within the subarachnoid and ventricular spaces. We modeled this with a sinusoidal functional based on [32]. With the oscillating function that is net positive, the gradient initially mimics this sinusoidal pattern but quickly converges to a gradient that averages ~13%.

Background concentration can negate bound CXCL12 gradient
As the single tumor cell is not the only contributor to CXCL12 gradients, it is useful to look at how background concentration can affect the development of the autologous gradient. After letting the cell establish its own gradient, we introduce CXCL12 to the background at 150 cs ( Figure 4A). The abrogation of the pericellular gradient can be seen over time as the background concentration completely overtakes any gradient development. This abrogation depends on a number of factors such as the bulk concentration of CXCL12 being released by the cell (4B), the background concentration of CXCL12 that the cell is surrounded by (4C), and the distance the cell is from the source of the concentration (4D). If the background CXCL12 is kept constant at 100 nM and the amount of CXCL12 the cell secretes is changed, there is a tipping point where the pericellular gradient starts to form -after 250 nM CXCL12 bulk concentration -after which it reaches a plateau (4B). If the reverse scenario is explored (constant cell secretion with changing background CXCL12), the pericellular gradient formation is abrogated around 1 mol/m 3 (4C). The position of the cell is also important where background concentration is concerned. As the cell is modeled farther from the source of the background concentration, the pericellular gradient becomes less disrupted to the point that at 70 µm from the source the autologous gradient is still observed for the length of time that the background concentration is applied (4D). If the cell is closer than 70 μm the background concentration completely abrogates pericellular gradient formation and actually causes a higher concentration of bound CXCL12 to develop upstream of the cell. If the background concentration is removed, the cell reaches a new set point of gradient formation which is notably lower than the starting %concentration the closer the cell is to the source.

The invading cell needs to be a certain distance from tumor border for gradient to develop
Looking into a more physiologically relevant parallel, cells that invade past the tumor border are subjected to varying CXCL12 levels based on the cells at the tumor border that may also secrete CXCL12 as shown in histological samples ( Figure 5A). By creating an approximate model with multiple cells on one side of the invading cell and applying CXCL12 secretion to a percentage of those cells, we can model the % concentration that an invading cell would have (5B, C). As the cell gets farther from the tumor border, the % concentration of bound CXCL12 around the cell increases. In other words, the pericellular gradient can develop without the interference of background CXCL12 secreted by the tumor border in areas where flow is outward. This is further impacted by the magnitude of fluid flow that the cell is under. At higher fluid velocities, the cell is actually able to develop a bound CXCL12 gradient closer to the tumor border. At the lower velocities, the cell cannot overcome the background CXCL12 even at 120 μm away from the border.

Cell type and morphology affects gradient formation
Up to this point, we have not utilized the full capability of the 3D model as we have concentrated on a readily understandable 2D outcome, %concentration. However, we can probe more complex parameters such as cell size and morphology. As shown in Supplemental Figure S7, cell shape and orientation has a large impact on gradient formation. As the cell increases in diameter, %concentration actually decreases. If we then imagine the cell as an ellipsoid, such as might be seen in a mesenchymal state, we can look into how orientation plays a role in gradient formation. We performed a rotational parametric sweep which shows %concentration is greatest along the long axis of the ellipsoid or when the cell is oriented in the direction of flow (Supplemental Figure S7). When introducing more complex geometry such as that of actual cells, gradient formation becomes heterogeneous and dependent on topography of the cell ( Figure 6). In order to accomplish this, we imaged different glioma cells that were within a 3D hydrogel (collagen and hyaluronic acid mixture) at 60x (S7A). The images were then rendered in 3D using a combination of ImageJ and meshing software (refer to the supplemental materials for a detailed description) and imported into the COMSOL model (S7B, C). The overall gradient stays largely the same at the pericellular level but the surface concentration that the cell sees is much more variable (S7D). This topographical heterogeneity is highly dependent on individual cell morphology which complicates matters when trying to understand gradient formation. It may well be that the overall concentration profile is more important than discrete locations on the cell, but this should not be dismissed out of hand.

Discussion
In this model, the circumstances leading to and causing autologous gradient formation are investigated through a variety of means in order to better understand forces potentially driving glioma invasion. These glioma cells are introduced to a variety of forces and tissue properties which are heterogeneous by nature. These can be from the environment such as stiffness of tissue, permeability and porosity of extracellular matrix, flow velocity magnitudes, and shear stress, or from internal origin such as regulation of cytokines, cytoskeletal remodeling, morphologic/phenotypic changes (epithelial to mesenchymal transition), secretion rates, or signaling kinetics. In Figure 1, we modulate individual parameters and observe the effects they have on the development of the chemokine gradient around the cell. Interestingly, velocity magnitudes seem to have the highest impact on this gradient formation. This may be unsurprising as one would expect higher flow rates to carry chemokines farther and there have been other published models which observe the same for similar chemokines and pathological states [25], [26], but it is also complicated by the binding kinetics and tissue properties that are inherent to the system. This does further strengthen the idea, though, that increased interstitial flow magnitudes like those seen in brain cancer can mediate an autologous chemotaxis response, adding to the invasion potential of tumor cells. The reaction kinetic coefficients were also investigated in order to get an idea of the effect that they have on the model, especially because there is a lack of experimental measurements of these values in the literature. While these should theoretically not change substantially (and not to the extent that we have manipulated in the model), it is worth looking at their impact. From the model results, the % concentration changes from -9 to 25%. This shows a large sensitivity in our model and one that should be further explored as more data is collected on the binding rates and reaction kinetics of chemokines and their various binding sites. Lastly, the diffusion coefficient is the other variable that the model is sensitive to, with changes to % concentration of -2 to 21%.
We have observed heterogeneous flow magnitude and direction in and around glioblastoma tumors [9]. In particular, there seems to be a role between flow direction and glioma progression. The MRI overlay on the histological section shown in Figure 2 depicts multidirectional flow on invading cells. The effects of these different flow vectors have yet to be explored on autologous chemotactic gradients. We show here that multidirectional flow has a direct impact on gradient formation, with the gradient aligning in the direction of the resultant flow profile. In instances where a cell might be seeing flow from opposite directions, the gradient remains mostly unaffected, but there will be other factors to consider such as how the resulting increase in pressure from this flow profile will impact the affected cell as there are other pathways by which flow could mediate tumor migration such as through mechanotransduction.
As mentioned previously, the brain has dynamic processes taking place all the time especially when considering fluid flow. It has been shown that fluid flow can fluctuate based on changes in respiration, heart rate, circadian rhythm, and body position (supine vs sitting) [33]- [35]. These flows happen at different time scales and rates, and are made even more complex when considering pathologic states such as cancer. There are many ways that flow can change within the context of cancer -directly due to pressure from the tumor and surrounding modified tissue, from drugs/drug delivery methods, or from surgical procedures (resection). Thus, it is useful to look at how gradient formation takes place over time and under these different circumstances to understand when/if autologous chemotaxis is a feasible method for cellular migration and the impact this might have on treatment outcomes. In the case of tumorigenesis there is an increase in pressure as the tumor grows, with heterogeneous flow magnitudes appearing. This can lead to a resulting scenario such as the ramped flow function shown in Figure 3. In the case of pulsatile flow, certain treatments such as CED cause fluctuations to interstitial fluid flow based on the treatment schema followed. The oscillating function mimics conditions that may be seen as a result of changing blood pressure or intracranial pressure from a number of things. From our model it seems that, given enough time, gradient formation eventually converges with small fluctuations at earlier time points.
Up to this point, we have only considered gradient formation in the context of a single cell secreting all of the chemokine. In reality, chemokines are released throughout tissue from different cell types within the microenvironment which could impact gradient formation if close enough to the tumor cell. In particular, this could come from tumor cells located at the tumor border or from cells that have been polarized by the tumor itself. If we let the model cell develop its gradient and then apply a background concentration of CXCL12, we see that the gradient around the cell is negated. There are multiple factors that impact this relationship between the background concentration and gradient formation, one of which is how much CXCL12 the cell is releasing. At some point, if the cell is secreting enough CXCL12 in relation to the background concentration, the gradient will still develop as shown in Figure 4. Depending on how close the cell is to the source of this background concentration, the gradient can actually form in the opposite direction (upstream of the cell instead of downstream), denoted by a negative % concentration value. If the cell is far enough away, the % concentration will be lessened, but actually seems to reform its gradient more slowly as seen by the slope of the line in Figure 4D.
If we introduce cells into the model to form a 'tumor border' and then place an 'invading cell' and modulate its distance from this border, we can see a similar trend where the gradient is abrogated when the cell is close to the border and starts to form at a specific distance from the border. The flow magnitude in the local environment also impacts this greatly, with increased velocity causing gradient formation to happen even when the cell is close to the tumor border. In the context of glioma, this is very interesting as it suggests several things. The first is that invading cells in areas of low flow magnitude must be using multiple modalities by which to move when first invading from the tumor border as the autologous gradient would not have developed in most cases. It is only once the cells reach a certain distance that autologous chemotaxis might be taking place to cause tumor cell movement. The other is that under higher flow magnitudes, cells could possibly undergo autologous chemotaxis even when near the tumor border as a way to invade into surrounding tissue. This seems more likely for flow magnitudes between 10 and 100 µm/s.
We also see through our model that cell shape, orientation, and sizes play a large role in gradient formation. These can all vary throughout the life of a cell and are certainly different between similar cells, not to mention between cell types. The model suggests that larger cells (25-50 µm) are less able to form gradients that would lead to migration due to autologous chemotaxis. This is of course when looking at a perfectly spherical cell. In instances where the cell is elongated as might occur in cells undergoing epithelial to mesenchymal transition or in cells that are activated or just possess such an inherent morphology, our model shows an increase in the % concentration when the cell is aligned in the direction of flow (its long axis is parallel with flow). This might have major implications for differences in cell migration based on phenotype/morphology (in terms of autologous chemotaxis), and also mean that some cells might be predisposed to autologous chemotaxis-effected movement. In order to probe deeper into the role of morphology on gradient formation, we introduced 3D-rendered cells (developed from confocal imaging) into the model. The concentration gradient changes based on topographical region of the cell, but looking at the pericellular gradient as a whole the heterogeneity is less pronounced. This could have consequences on autologous chemotaxis based on which factor contributes more to the chemotactic effect, but is not in the scope of this paper. This also points to a heterogeneity between cells of the same type based on morphology as well as between cell types. More cells will need to be analyzed in order to see any apparent trends.
There are many future directions and potentials for this model. The obvious first question is what amount of CXCL12 and what gradient is needed to cause chemotaxis to occur in glioma cells? In order to understand this, in vitro studies will need to be carried out that vary the concentration of CXCL12 across CXCR4 positive cells. More info is also needed on the specifics of the protease release reaction kinetics and the concentration of HS in the ECM. Having these in hand will drastically increase the robustness of the model and allow us to more confidently interpret the model results. In future, the model could also be improved by introduction of fluid-structure interactions to account for other fluidmediated effects, such as shear stress and cell morphology changes, on the formation of chemokine gradients.

Conclusions
In conclusion, we have explored the impact of physiological flows and parameters on the formation of autologous gradients around tumor cells. We see that in every case, gradients can form, however, these are often reduced as more tissue-level relevance is added. Not only is the nature of the flow important but also the extracellular context including background gradients and other cells as would be present at the tumor border.
Further, the cell itself, when modeled based on an actual cell which is non-spherical in shape, we see differential gradients depending upon the 2-dimensional cross-section which is examined. These explorations tell us that the formation of autologous gradients in vivo is much more complex than what has been modeled before and what has been replicated in vitro.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1: Schematic of computational model baseline. Figure S2: Mesh refinement analysis Figure S3: Sensitivity analysis of krel Figure S4: Parametric sweeps of input variables and their corresponding effects on %concentration. Figure S5: Additional multidirectional flow modeling. Figure S6: Time effect on gradient formation for baseline condition. Figure S7: Cell size and orientation impact gradient formation. Supplemental Methods.  Data Availability Statement: Data is available upon request by contacting the corresponding author.