Computational Simulation of Exosome Transport in Tumor Microenvironment

Cellular exosome-mediated crosstalk in tumor microenvironment (TME) is a critical component of anti-tumor immune responses. In addition to particle size, exosome transport and uptake by target cells is influenced by physical and physiological factors, including interstitial fluid pressure, and exosome concentration. These variables differ under both normal and pathological conditions, including cancer. The transport of exosomes in TME is governed by interstitial flow and diffusion. Based on these determinants, mathematical models were adapted to simulate the transport of exosomes in the TME with specified exosome release rates from the tumor cells. In this study, the significance of spatial relationship in exosome-mediated intercellular communication was established by treating their movement in the TME as a continuum using a transport equation, with advection due to interstitial flow and diffusion due to concentration gradients. To quantify the rate of release of exosomes by biomechanical forces acting on the tumor cells, we used a transwell platform with confluent triple negative breast cancer cells 4T1.2 seeded in BioFlex plates exposed to an oscillatory force. Exosome release rates were quantified from 4T1.2 cells seeded at the bottom of the well following the application of either no force or an oscillatory force, and these rates were used to model exosome transport in the transwell. The simulations predicted that a larger number of exosomes reached the membrane of the transwell for 4T1.2 cells exposed to the oscillatory force when compared to controls. Additionally, we simulated the interstitial fluid flow and exosome transport in a 2-dimensional TME with macrophages, T cells, and mixtures of these two populations at two different stages of a tumor growth. Computational simulations were carried out using the commercial computational simulation package, ANSYS/Fluent. The results of this study indicated higher exosome concentrations and larger interstitial fluid pressure at the later stages of the tumor growth. Quantifying the release of exosomes by cancer cells, their transport through the TME, and their concentration in TME will afford a deeper understanding of the mechanisms of these interactions and aid in deriving predictive models for therapeutic intervention.

Introduction

The tumor microenvironment (TME), which contains multiple cell types, blood and lymph vessels, and extracellular matrix (ECM), is often immunosuppressive, blocking anti-tumor immunity, and mediating drug resistance (1). Previous studies have shown that tumors have significantly increased hydrostatic pressure and solid force (tension, compression, and shear), and these forces activate signaling pathways and drive changes in tumor cell proliferation, adhesion, migration, and survival (2). Emerging evidence indicates that biomechanical forces in the TME may influence tumor progression due to anti-tumor immune response through cancer cell-immune cell crosstalk, including that mediated by vesicular transport (3, 4). As a key mediator of cell-cell signaling within the TME, tumor cell exosomes have emerged as important regulators of several aspects of tumor progression (5), including tumor initiation, progression and invasion, and preparation of pre-metastatic niches (6, 7). Exosomes are 50–150 nm membrane-bound extracellular vesicles (EVs) released by various cells in the TME, including tumor cells (8). Exosomes contain proteins and microRNAs that can impact immune cell function. We and others have demonstrated that tumor-derived exosomes promote macrophages polarization (9–11). They contribute to immune regulation and can create an immuno-privileged environment within tumors (8, 12).

Exosomes secreted by the primary tumor lesion have been shown to aid in the formation of metastatic lesion in distant tissues (i.e., pre-metastatic niche formation) before the cancer cells themselves migrate to the particular tissue. In 2012, Peinado et al. (13) showed that when exosomes derived from melanoma cells, with highly metastatic potential to the lungs, were administered to mice, they promoted vascular permeability through the introduction of bone marrow cells to the lungs and thereby contributed to the formation of a pre-metastatic niche. Similarly, pancreatic cancer cell-derived exosomes formed a pre-metastatic niche for liver metastasis (6, 14). Internalization of pancreatic cancer cell-derived exosomes containing macrophage migration inhibitory factor by hepatic Kuppfer cells, induced fibronectin and TGF-β production by hepatic stellate cells and ultimately formed a pre-metastatic niche through introduction of bone marrow cells to the liver (15). In the context of infections, exosome transfer by lymphatic flow from the periphery to the lymph node has been proposed to be a mechanism for rapid exchange of infection-specific information that precedes the arrival of migrating cells, thus priming the node for a more effective immune response (16). Recent studies have also investigated the utility of exosomes for drug delivery and as carriers of bioactive molecules (17, 18). In vivo studies demonstrate that after being internalized at the administration site and transferred to the systemic circulation, exosomes pass blood–tissue barriers and arrive in each tissue (17–19). Intravenously administered exosomes have been shown to be mainly distributed in the organs with a mononuclear phagocyte system such as the liver, spleen, lungs, and kidneys (6, 13, 20). In our own studies, intratumoral injection of tumor-derived exosomes resulted in internalization by mononuclear phagocytes including immune suppressive myeloid-derived suppressor cells (21). The release of exosomes by cancer cells, their transport through the TME, and their crosstalk with immune cells once they reach the tissue are not well-understood. Even though the importance of exosomes in TME is well-documented in literature, tools to quantify exosomes in TME are lacking. This necessitates the development of a numerical approach to simulate the transport of exosomes and estimate concentration gradients of exosomes in TME.

With the advent of faster and better computational resources and algorithms, researchers have modeled various aspects of the TME including simulation of tumor growth in the microenvironment (22–25). Several studies have comprehensively reviewed the existing models for tumor growth, tumor progression and morphology (22, 23). These models took into account the genetic characteristics of the tumor and the TME. Bresch et al. (24) used partial differential equations to model tumor growth and to estimate tumor densities, and level set methods to model the membrane. Sciumè et al. (25) developed an approach to model the TME with tumor and healthy cells, interstitial fluid, and extracellular matrix using continuum mechanics principles. Crespo et al. documented the importance of mathematical modeling of heterogeneous systems in tumor microenvironments and interplay of these elements in TME (26). Further, they elucidated the use of computational simulations in personalized cancer therapy and precision medicine.

Various studies have addressed the computational simulation of interstitial fluid flow and estimation of shear stress and fluid pressure (27–30). Mitchell and King (27) recognized that mechanical forces such as fluid shear stress can influence cancer metastasis. They used the Stokes equation, Darcy's law, and Brinkman equation to estimate the fluid properties in TME. Welter and Rieger (29) postulated that the elevated interstitial fluid pressure due to a tumor could influence the interstitial fluid flow and delivery of drugs and nutrients to the cells. They modeled TME as a porous medium and the movement of chemicals in TME was modeled using a transport equation. A non-invasive technique to estimate the response of cancer treatment using imaging, pharmacokinetics, and interstitial fluid flow modeling was presented by Swinburne et al. (30). In this modeling, the governing equations were taken as the Navier-Stokes equation and the extracellular matrix was modeled as a porous medium. Voronov et al. (31) used a numerical approach to estimate the fluid shear stress in a scaffold. A high resolution micro-CT was used to extract the detailed geometry of the scaffold and the lattice Boltzmann method used for the modeling of fluid flow inside the scaffold. Kim et al. used a medical imaging technique (32) to extract high fidelity 3-D geometry of a tumor and mathematical models to delineate blood flow and molecular transport. The goal of this multi-scale modeling was to develop a predictive model for tumor angiogenesis. Another approach to model TME are cell or agent-based methods, which simulate the dynamic evolution of individual entities called agents based on a set of rules (33–36). A survey of different cell-based mathematical models to analyze TME and their relative strengths and weaknesses are presented in Rejniak and McCawley (37). Cell-based approaches are particularly useful in modeling individual cell behavior such as polarization of immune cells, which provide insight into the dynamic processes in heterogeneous TME.

Although various aspects of TME have been modeled using computational simulations, the computational tools to simulate exosome transport in TME are lacking. Quantifying exosome concentration and factors that regulate their interactions with various cells in the TME will afford a deeper understanding of the mechanisms of interaction and aid in deriving predictive models for therapeutic intervention. This article focuses on a computational technique to model the transport of exosomes in the TME by taking into account exosome production by cancer cells, exosome advection due to interstitial flow, and exosome diffusion due to concentration gradients.

Materials and Methods

Computational Modeling

In this study, because of the sub-micron size of exosomes, their movement in the TME was treated as a continuum and modeled using an advection-diffusion transport equation, with advection due to interstitial flow and diffusion due to concentration gradients (38). In general, an exosome transport equation can be written as,

$\begin{array}{ll}\frac{\partial c_k}{\partial t}+\frac{\partial u_i{c}_k}{\partial x_i}=\frac{\partial }{\partial x_i}\left({\Gamma }_k\frac{\partial c_k}{\partial x_i}\right)+S_{c_k},k=1,\dots ,N& (1)\end{array}$

Where c_k is the concentration of the k^th type of exosome in a heterogeneous field with N number of exosome types, u_i is the interstitial fluid velocity components, Γ_k is the diffusion coefficient, and S_ck is the exosome source or sink term. In the above equation, the second term on the left-hand side models the movement of exosomes due to interstitial fluid flow, the first term on the right-hand side models the diffusion of exosomes due to concentration gradients, and the last term on the right-hand side represents the production of exosomes by donor cells (tumor cells) and uptake of exosomes by receptor cells (immune cells). In this study, the uptake of exosomes by immune cells were not taken into consideration. The diffusion coefficient for the simulation was taken from the data published in the literature (38–42). The interstitial fluid flow in TME is governed by the Navier-Stokes equations (43–45), which is the mathematical representation of the conservation of mass, momentum, and energy. The fluid flow in TME is incompressible in nature, and heat transfer is also not an important factor for the flow in TME (46). Therefore, the energy equation was not considered for the following simulations. The turbulent effects were not modeled due to the fact that the Reynolds number for the interstitial flow was very small. Therefore, the laminar flow assumption was used for the simulations.

Computational simulations presented in this article were carried out using the commercial computational simulation package called ANSYS/Fluent (47). It has capabilities to model Navier-Stokes and transport equations for different flow regimes using various numerical schemes, and is commonly used in academia and industries. In the present simulations, the transport of exosome in TME was carried out in two stages using a staggered approach. In the first stage the Navier-Stokes equations were solved to estimate the velocity field. In the second stage, the transport equations were solved using the velocity field predicted from the first step. In these simulations, it was assumed that the transport of exosomes in TME would not affect the underlying velocity field. Geometries for TME, presented in this article were extracted from the images created using Biorender (48). These images were segmented using edge detection techniques to get the outer boundaries of tumor and immune cells present in TME, with the application of the geometry modeling program SpaceClaim available in Ansys (47).

Validation of the Numerical Approach

Validation of the numerical methods using experimental data or analytical solution of an appropriate problem is an important step in computational simulations. The numerical approach to solve the advection-diffusion equation used in this study was validated using a Gaussian source in a planar shear flow (49, 50). The computational domain for this simulation was taken as a square with dimension 100 × 100 mm, with the origin at the geometric center. In this test case, the y-component of the velocity was taken as zero and the x-component of the velocity was assumed to vary linearly as U + λy, where U and λ are constants. At time equals zero, the concentration of the transport variable was taken as an impulse function at the origin, multiplied by a constant M. The exact solution of the concentration for this Gaussian source in a planar shear flow at any spatial location and at any time is given by the relation (49, 50),

$\begin{array}{l}c(x,y,t)=\\ \frac{M}{4\pi Dt{\left(1+\frac{{(\lambda t)}^2}{12}\right)}^{\frac{1}{2}}}exp\left\{-\frac{1}{4Dt}\left[y^2+\frac{{(x-\left(U-0.5\lambda y\right)t)}^2}{\left(1+\frac{{(\lambda t)}^2}{12}\right)}\right]\right\}& (2)\end{array}$

Where c is the concentration and D is the diffusion coefficient. More details on the analytical solution of this benchmark test case can be found in references (42, 43). Numerical values of different parameters used in this benchmark test case are listed in Table 1. To keep the concentration as finite values in the computational domain at the beginning of the simulation, the exact solution at time equals 1 min was used.

TABLE 1

Table 1. Parameters used for the benchmark test case.

Cell Culture and Exposure to Strain

The rate of release of exosomes is influenced by biomechanical forces acting on the tumor cells, and to quantify this, an experimental setup depicted in Figure 1 was used. The murine TNBC cell line 4T1.2 (an aggressive clone derived from 4T1) was obtained from Dr. Robin L. Anderson's laboratory (Peter McCallum Cancer Center, Australia). 4T1.2 cells were cultured in Dulbecco's Modified Eagle Medium (DMEM) supplemented with 10% FBS and 10 mM HEPES (MP Biomedicals, Santa Ana, CA). Prior to the exposure to tensile strain, the 4T1.2 cells were stained with the lypophilic dye PKH26 Red Fluorescent Cell Linker (Sigma, St. Louis, MO), per the manufacturer's instructions. 2.5 × 10⁵ 4T1.2 were seeded on collagen coated 6 well UniFlex culture plates (Flexcell International Corporation, Burlington, NC) and cultured to confluence. Once confluent, the media was changed to exosome depleted growth media and the plates were subjected to 10% uniaxial oscillatory strain at 0.3 Hz for 48 h, 10% constant strain for 48 h, or no strain for 48 h using a FlexCell FX-6000 or FX-5000 Tension System.

FIGURE 1

Figure 1. Schematic of experimental protocol. (A) 4T1.2 cells were seeded in BioFlex plates (1), once cells reached confluence they were exposed to OF strain (or NF) (2). Strain exposed (or NF) cells were either put into co-cultured with immune cells (T cells and macrophages) using transwell (3) or conditioned media was collected for exosome isolation. (B) Exosome release rate quantified following OF or NF. n = 8 per condition from 3 experiments, 2–3 replicates per experiment. (C) Exosome internalization by T cells and macrophages was quantified using Imagestream flow cytometry, n = 2–3 per condition.

Preparation of Exosome-Depleted Media

The exosome-depleted medium was prepared as previously described (51). Briefly, DMEM media supplemented with 20% FBS was centrifuged using an ultracentrifuge overnight at 100,000 × g at 4°C. The supernatant was filtered through a 0.2 μm cellulose acetate filter (Corning, NY). The exosome depleted media was then diluted 1:1 with DMEM to make a final concentration of 10% FBS and 10 mM HEPES.

Isolation of Exosomes From Conditioned Media and Quantification of Concentration

Extracellular vesicles (EVs) were isolated as previously described (21). Briefly, conditioned media was collected following 48 h exposure to NF or OF and centrifuge at 2,000 × g to remove any cell fragments or apoptotic bodies. The supernatant was then incubated with the Total Exosome Isolation Reagent for Cell Culture Media kit (ThermoFisher, Waltham, MA) per the manufacturer's protocol. Purified exosomes were stored in 50 μl of PBS at −80°C. Exosome concentration was quantified using Imagestream flow cytometry.

Immune Cell Exosome Internalization

0.4 μm pore size transwell filters (MilliporeSigma, Burlington, MA) were added on top of strained or control 6 well plates containing PKH26⁺ 4T1.2 cells. Subsequently, 5 × 10⁴ purified naïve lung macrophages (CD11b⁺CD64⁺), or 3 × 10⁵ purified splenic T cells (CD45⁺CD8⁺), or a mixture of these cell populations (1:3 ratio of macrophages to T cells, total of 2 × 10⁵) were added into the top well of the transwell filter. After 24 h, cells on the upper surface of the transwell filter were removed and tumor cell exosome internalization was evaluated using ImageStream flow cytometry, evaluating cell specific markers CD8 and CD11b, along with PKH26, to indicate exosome internalization.

Results

Results From the Validation of the Numerical Approach

For the simulation of the validation test case, the concentration of the transport variable at the beginning of the simulation is shown in Figure 2A, and the concentration after 45 min is shown in Figure 2B. These data show the advection due to the shear flow and diffusion due to the concentration gradient. Also, it can be seen from the result that as time progresses, the concentration gets distorted due to the velocity gradient normal to the flow direction. Computations were carried out to calculate the concentrations at 45 min using different numerical schemes and a quantitative comparison of the computed results with the analytical data is presented in Figure 3.

FIGURE 2

Figure 2. The initial conditions used for the benchmark simulations and the predicted concentration at the end of the simulation. In these figures, red represents higher concentrations and blue represents lower concentrations. (A) Initial concentration used at the beginning of the simulation. (B) Computed concentration after 45 min.

FIGURE 3

Figure 3. Comparison of computed concentrations using different numerical schemes with the analytical solutions. Simulations were conducted using first order, second order upwind, second order QUICK, and third order MUSCL schemes for the discretization of the spatial derivatives in the transport equation (37–39). Computed concentrations are compared with the analytical results along the horizontal and vertical lines shown in Figure 2B. (A) Comparison of the concentration along the horizontal line. (B) Zoomed in view of the comparison in the rectangle marked in (A). (C) Comparison of the concentration along the vertical line. (D) Zoomed in view of the concentration in the rectangle marked in (C).

3-D Transwell Modeling

To mimic the experimental setup depicted in Figure 1, a 3-D computational model of the transwell geometry was created, as shown in Figure 4A, and was used for the simulation of exosome transport. This geometry was created using the computer aided geometry (CAD) modeler available in ANSYS (47). The cells that release exosomes were placed at the bottom chamber and the release rates of the exosomes from these cells were taken from the experiments described in the previous section. The experimental measurements described in section Isolation of Exosomes from Conditioned Media and Quantification of Concentration showed that the average exosome release rates were 515.3, and 986.3 per hour per 10⁵ cells, for NF and OF groups, respectively (Figure 1B). In co-culture experiments, internalization of exosomes by immune cells could reduce the exosome concentration near these cells, which could potentially affect the diffusion. However, experimental evidence showed that the concentration change due to internalization was small (Figure 1C) and therefore, exosome internalization measurements from the experiment were assumed negligible and was not modeled. The presented numerical approach will be further refined by adding a sink term in the transport equation (equation 1) to represent exosome internalization in future studies. For the computational simulations, 3.6 × 10⁶ cells were assumed to be placed at the bottom of the transwell for both NF and OF groups, and the exosome source fluxes from the bottom surface were estimated using the exosome release rates from the experiments and the area of the bottom surface of the transwell.

FIGURE 4

Figure 4. Computational domain and the results from modeling of exosome transport in a transwell. For this modeling, exosome release rate for the cells placed at the bottom surface is taken as 515.3 and 986.3/h./10⁵ cells for NF and OF, respectively. (A) A cross-sectional view of the 3-dimensional computational domain used for the simulation. (B) Exosome gradient on a plane through the middle of the transwell for NF after 72 h. (C) A comparison of the number of exosomes at the membrane for NF and OF groups as a function of time.

For these simulations, the fluid medium in the transwell arrangement was taken as stagnant. This resulted in the transport of exosomes in the transwell based completely on diffusion. Therefore, the velocity field was set as zero and the Navier-Stokes equations were not used to calculate the flow field. The diffusion coefficient for this simulation was taken as 2.0 × 10⁻³ m²/s. The unsteady transport equation was used to estimate the exosome concentration in the computational domain. The initial concentration of the exosomes in the transwell was taken as zero and the computational simulation calculated the time evolution of exosome concentration in the transwell. Simulations were carried out for 72 h and the concentration of exosomes in the computational domain for NF and the comparison of the number of exosomes that reached the membrane for both NF and OF groups are shown in Figures 4B,C.

Modeling the Tumor Microenvironment

The presence of tumor in a microenvironment not only influences the dynamic behavior of the cells in the microenvironment, but also changes the interstitial fluid flow properties such as pressure and shear stress distributions, and flow patterns. As the tumor grows, obstruction to flow increases and the exosome release rate also increases due to the larger number of cancer cells. To model TME, two different configurations at two different stages of the tumor growth in TMEs with macrophages, T cells, and MDSCs were created using Biorender (48). This approach of generating geometries for TME using Biorender images was considered for this proof-of-principle study. However, more realistic TME from medical images will be considered in follow-up studies. The TME with different cell types at the later stage of growth is shown in Figure 5A. The distribution and proportions of these cell populations were based on what is observed normally in the TME in murine models of both lung and breast cancer (21, 52–55). The extracted geometry using image segmentation techniques and a zoomed-in view of the computational mesh near the tumor surface are shown in Figures 5B,C, respectively. In this simulation, the interstitial fluid was assumed to be coming from the left side of the computational domain with a velocity of 0.75 μm/s. The Navier-Stokes equations together with laminar flow assumptions were used for the calculation of the velocity field in TME. This computed velocity field was used for the evaluation of the advective component of the transport equation. The exosome release rates for the initial and final stages of the tumor were taken as 2.5 and 7.5 exosomes per second, respectively. The computed streamlines in TME and the concentration of exosomes after 45 min are shown in Figure 6.

FIGURE 5

Figure 5. Biorender image, computational domain, and a view of the mesh for the tumor microenvironment with macrophages, T cells, and MDSCs. The tumor microenvironments used for the simulations are created using Biorendor (48), and the rendered images is segmented to extract the geometry of various cell types in the computational domain. (A) Rendered image using Biorender. (B) Extracted computational domain using image processing techniques. (C) A zoomed in view of the mesh in TME used for the computational simulation.

FIGURE 6

Figure 6. Results from the simulation of interstitial fluid flow and exosome concentration in TME at two stages in tumor growth. Exosomes produced by the tumor cells at a constant rate are carried by the interstitial fluid flow and diffused based on the concentration gradient. (A) Predicted exosome concentration in an early stage. (B) Predicted exosome concentration in a later stage. The predicted concentration of exosomes is ~3.5 times greater for the later stage as compared to the early stage.

Discussion

The validation of the numerical approach for the computational modeling was carried out using the advection and diffusion of a material in a planar shear flow, which has an analytical solution. The order of the discretization of the spatial derivative terms in the transport equation (equation 1) influences the accuracy of computational simulations. To estimate the accuracy of these discretizations and to select an appropriate numerical scheme for the rest of the simulations presented in this paper, simulations were conducted using first order, second order upwind, second order QUICK, and third order MUSCL schemes (43–45). The computed concentrations at 45 min from these simulations were compared with the exact solution (equation 2) in Figure 3. This comparison shows that the results from the second and third order schemes are in very good agreement with the analytical results, and the selected schemes are able to model both advection and diffusion accurately. Therefore, for the rest of the calculations, a second order numerical scheme was used.

The transwell simulations modeled the exosome transport through a stagnant region due to the concentration gradient. As expected, exosome gradient on a plane showed a higher concentration of exosomes at the bottom wall and a lower concentration at the top of both bottom and top wells. In this figure, red represents a higher concentration of exosomes and blue represents a lower concentration. Figure 4C compares the number of exosomes that have reached the membrane at the base of the top well for both NF and OF groups. Computational simulations showed that at all times, a larger number of exosomes reached the membrane for the OF group compared to the NF group. After 72 h, ~2.86 × 10⁴ exosomes reached the membrane for 4T1.2 cells exposed to OF compared 1.58 × 10⁴ exosomes for those exposed to NF. The difference in the number of exosomes passing through the membrane is due to the higher release rate for the OF group, resulting in a higher concentration gradient and faster diffusion. The rate of release of exosomes by the tumor cells in the OF group is ~1.9 times compared to the release rate for the ones in the NF group. This shows the number of exosomes that reached the membrane is close to being linearly proportional to the exosome release rate.

As shown in Figure 6, the release rate for the later stage of the tumor was taken as a larger value due to the presence of a larger number of cells in the tumor. Here, red represents higher and blue represents lower exosome concentrations. In the later stage (B), the predicted average concentration of the exosomes was ~3.5 times greater than the predicted exosome concentration in the early stage (A). The streamline patterns in the TME for the two different stages of the tumor growth is also plotted in Figure 6. It can be seen from the figure that the flow pattern in the TME is more complex for the later stage of the tumor growth, and this is due to more blockage for the interstitial fluid flow. In addition to the complex flow pattern, the interstitial fluid pressure was higher for the later stages of the tumor growth compared to the early stage, see Figure 7. This could also lead to a higher rate of exosome release due to the larger mechanical forces acting on the tumor surface. However, in this simulation the influence of higher interstitial fluid pressure on exosome release was not taken into account.

FIGURE 7

Figure 7. Pressure distribution in the TME at the early and later stages of tumor growth. The Navier-Stokes equations are used to estimate the interstitial fluid flow velocities and the pressure distribution within TME. (A) Predicted pressure distribution at early stage. (B) Predicted pressure distribution at later stage. The predicted pressure distribution is greater for the later stage of the tumor growth, due to greater obstruction for the flow by a larger tumor.

The concentration of exosomes in TME is an important parameter in cellular crosstalk and evolution of the dynamic environment. Drug delivery and microRNA delivery via exosomes for specific targets in the TME is an active area of investigation that has significant translational potential. This paper presents a numerical approach to model the transport of exosomes in the TME and quantify exosome concentrations. Results are presented to validate the approach using a benchmark test case, and the accuracy of different numerical schemes are presented by comparing numerical results with the exact solution. The numerical approach presented in our studies is a valuable tool to quantify exosome concentration, exosome gradient, and time evolution of exosome concentration in a TME. This approach can be combined with agent-based models to simulate exosome uptake by immune cells and their polarization. This will enable the evaluation of the influence of different parameters, including the magnitude of force, the frequency and duration of the application, and distance from the tumor cells on the impact of exosomes on immune populations and their polarization rather quickly, without conducting a large number of experiments. This will also make it easier, faster, and more cost-effective to study the effect of the wide range of these parameter values on immune cell polarization that determine the effector or suppressor function at the TME. The use of this numerical approach for in-vivo applications requires the proper estimation of tissue properties in a heterogeneous environment, accurate representation of distribution of various cell types in the dynamic ECM environment, and the uptake and polarization of immune cells. However, this approach could serve as a tool to study the propagation of exosomes secreted by the primary tumor lesion and pre-metastatic niche formation in distant tissues, before the cancer cells themselves migrate to the particular tissue, and this will be explored in future studies. Additionally, modeling the dynamics of exosome transport in in-vitro models will provide insights to enhance the understanding of exosome transport in vivo from the primary tumor site to metastatic sites of the tumor and how exosome transport may influence metastatic niche.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Ethics Statement

All cell lines used in the study were screened for mycoplasma and cultured following recommended guidelines.

Author Contributions

RK, BT, and JB designed and conducted the simulations. KG, YW, SP, JB, and JD designed and conducted the experiments. KG and JB served as the interface to transfer information to and fro between computations and experiments. All authors contributed toward preparation of the manuscript and revisions.

Funding

These studies were supported by Breast Cancer Research Foundation of Alabama awarded to JB and JD. A part of this study was sponsored by the National Institute of Health under the grant R01 CA 184770 to SP.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The authors would like to acknowledge the support from UAB Research Computing for the facilities and the software used in this research.

Abbreviations

CAD, Computer aided geometry; DMEM, Dulbecco's Modified Eagle Medium; ECM: extracellular matrix; EVs, Extracellular vesicles; MDSCs, Myeloid-derive suppressor cells; MUSCL, Monotonic Upstream-centered Scheme for Conservation Laws; NF, No force; OF, Oscillatory force; QUICK, Quadratic Upstream Interpolation for Convective Kinematics; TME, Tumor microenvironment; TNBC, Triple Negative Breast Cancer.

References

1. Hui L, Chen Y. Tumor microenvironment: sanctuary of the devil. Cancer Lett. (2015) 368:7–13. doi: 10.1016/j.canlet.2015.07.039

CrossRef Full Text | Google Scholar

2. Jain RK, Martin, Stylianopoulos JD. The role of mechanical forces in tumor growth and therapy. Annu Rev Biomed Eng. (2014) 16:321–46. doi: 10.1146/annurev-bioeng-071813-105259

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Sottnik JL, Dai J, Zhang H, Campbell B, Keller ET. Tumor-induced pressure in the bone microenvironment causes osteocytes to promote the growth of prostate cancer bone metastases. Cancer Res. (2015) 75:2151–8. doi: 10.1158/0008-5472.CAN-14-2493

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Seager RJ, Hajal C, Spill F, Kamm RD, Zaman MH. Dynamic interplay between tumour, stroma, and immune system can drive or prevent tumour progression. Converg Sci Phys Oncol. (2017) 3:034002. doi: 10.1088/2057-1739/aa7e86

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hoshino A, Costa-Silva B, Shen TL, Rodrigues G, Hashimoto A, Tesic Mark M, et al. Tumour exosome integrins determine organotropic metastasis. Nature. (2015) 527:329–35. doi: 10.1038/nature15756

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Costa-Silva B, Aiello NM, Ocean AJ, Singh S, Zhang H, Thakur BK, et al. Pancreatic cancer exosomes initiate pre-metastatic niche formation in the liver. Nat Cell Biol. (2015) 17:816–26. doi: 10.1038/ncb3169

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Becker A, Thakur BK, Weiss JM, Kim HS, Peinado H, Lyden, et al. Extracellular vesicles in cancer: cell-to-cell mediators of metastasis. Cancer Cell. (2016) 30:836–48. doi: 10.1016/j.ccell.2016.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Kurywchak P, Tavormina J, Kalluri R. The emerging roles of exosomes in the modulation of immune responses in cancer. Genome Med. (2018) 10:23. doi: 10.1186/s13073-018-0535-4

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Pritchard A, Tousif S, Wang Y, Hough K, Khan S, Strenkowski J, et al. Lung tumor cell-derived exosomes promote M2 macrophage polarization. Cells. (2020) 9:5. doi: 10.3390/cells9051303

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Popena I, Abols A, Saulite L, Pleiko K, Zandberga E, Jekabsons K, et al. Effect of colorectal cancer-derived extracellular vesicles on the immunophenotype and cytokine secretion profile of monocytes and macrophages. Cell Commun Signal. (2018) 16:1. doi: 10.1186/s12964-018-0229-y

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Shinohara H, Kuranaga Y, Kumazaki M, Sugito N, Yoshikawa Y, Takai T, et al. Regulated polarization of tumor-associated macrophages by miR-145 via colorectal cancer-derived extracellular vesicles. J Immunol. (2017) 199:1505–15. doi: 10.4049/jimmunol.1700167

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Harris DA, Patel SH, Gucek M, Hendrix A, Westbroek W, Taraska W, et al. Exosomes released from breast cancer carcinomas stimulate cell movement. PLoS ONE. (2015) 10:e0117495. doi: 10.1371/journal.pone.0117495

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Peinado H, Aleckovic M, Lavotshkin S, Matei I, Costa-Silva B, Moreno-Bueno G, et al. Melanoma exosomes educate bone marrow progenitor cells toward a pro-metastatic phenotype through MET. Nat Med. (2012) 18:883–91. doi: 10.1038/nm.2753

CrossRef Full Text | Google Scholar

14. Yu Z, Zhao S, Ren L, Wang L, Chen Z, Hoffman RM, et al. Pancreatic cancer-derived exosomes promote tumor metastasis and liver pre-metastatic niche formation. Oncotarget. (2017) 8:63461–83. doi: 10.18632/oncotarget.18831

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Matsushita H, Yang YM, Pandol SJ, Seki E. Exosome migration inhibitory factor as a marker and therapeutic target for pancreatic cancer. Gastroenterology. (2016) 150:1033–5. doi: 10.1053/j.gastro.2016.02.051

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Srinivasan S, Vannberg FO, Dixon JB. Lymphatic transport of exosomes as a rapid route of information dissemination to the lymph node. Sci Rep. (2016) 6:24436. doi: 10.1038/srep24436

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Takahashi Y, Nishikawa M, Shinotsuka H, Matsui Y, Ohara S, Imai T, et al. Visualization and in vivo tracking of the exosomes of murine melanoma B16-BL6 cells in mice after intravenous injection. J Biotechnol. (2013) 165:77–84. doi: 10.1016/j.jbiotec.2013.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Kim MS, Haney MJ, Zhao Y, Yuan D, Deygen I, Klyachko NL, et al. Engineering macrophage-derived exosomes for targeted paclitaxel delivery to pulmonary metastases: in vitro and in vivo evaluations. Nanomedicine. (2018) 14:195–204. doi: 10.1016/j.nano.2017.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Das CK, Jena BC, Banerjee I, Das S, Parekh A, Bhutia SK, et al. Exosome as a novel shuttle for delivery of therapeutics across biological barriers. Mol Pharm. (2019) 16:24–40. doi: 10.1021/acs.molpharmaceut.8b00901

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yuan D, Zhao Y, Banks WA, Bullock KM, Haney M, Batrakova E, et al. Macrophage exosomes as natural nanocarriers for protein delivery to inflamed brain. Biomaterials. (2017) 142:1–12. doi: 10.1016/j.biomaterials.2017.07.011

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wang Y, Goliwas KF, Severino PE, Hough KP, Van Vessem D., Wang H, et al. Mechanical strain induces phenotypic changes in breast cancer cells and promotes immunosuppression in the tumor microenvironment. Lab Invest. (2020) 100:1503–16. doi: 10.1158/1538-7445.AM2020-1711

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Roose T, Chapman SJ, Maini P.K. Mathematical models of avascular tumor growth. SIAM Rev. (2007) 49:179–208. doi: 10.1137/S0036144504446291

CrossRef Full Text | Google Scholar

23. Macklin P, Lowengrub J. Nonlinear simulation of the effect of microenvironment on tumor growth. J Theor Biol. (2007) 245:677–704. doi: 10.1016/j.jtbi.2006.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Bresch D, Colin T, Grenier E, Ribba B, Saut O. Computational Modeling of solid tumor growth: the avascular stage. SIAM J Sci Comput. (2010) 32:2321–44. doi: 10.1137/070708895

CrossRef Full Text | Google Scholar

25. Sciumè G, Gray WG, Ferrari M, Decuzzi P, Schrefler BA. On computational modeling in tumor growth. Arch Comput Methods Eng. (2013) 20:327–52. doi: 10.1007/s11831-013-9090-8

CrossRef Full Text | Google Scholar

26. Crespo I, Coukos G, Doucey MA, Xenarios I. Modelling approaches to discovery in the tumor microenvironment. Cancer Immunol Ther. (2018) 1:26–40. doi: 10.35841/cancer-immunology.1.1.26-40

CrossRef Full Text | Google Scholar

27. Mitchell MJ, King, Computational MR. and experimental models of cancer cell response to fluid shear stress. Front Oncol. (2013) 3:44. doi: 10.3389/fonc.2013.00044

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kwon RY, Frangos JA. Quantification of lacunar-canalicular interstitial fluid flow through computational modeling of fluorescence recovery after photobleaching. Cell Mol Bioeng. (2010) 3:296–306. doi: 10.1007/s12195-010-0129-8

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Welter M, Rieger H. Interstitial fluid flow drug delivery in vascularized tumors: a computational model. PLoS ONE. (2013) 8:e70395. doi: 10.1371/journal.pone.0070395

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Swinburne N, LoCastro E, Paudyal R, Oh JH, Taunk NK, Shah A, et al. Computational modeling of interstitial fluid pressure and velocity in non-small cell lung cancer brain metastases treated with stereotactic radiosurgery. Front Neurol. (2020) 11:402. doi: 10.3389/fneur.2020.00402

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Voronov R, Vangordon S, Sikavitsas VI, Papavassiliou DV. Computational modeling of flow-induced shear stresses within 3D salt-leached porous scaffolds imaged via micro-CT. J Biomech. (2010) 43:1279–86. doi: 10.1016/j.jbiomech.2010.01.007

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Kim E, Stamatelos S, Cebulla J, Bhujwalla ZM, Popel AS, Pathak, et al. and computational modeling of blood flow in the tumor vasculature. Ann Biomed Eng. (2012) 40:2425–41. doi: 10.1007/s10439-012-0585-5

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Metzcar J, Wang Y, Heiland R, Macklin P. A review of cell-based computational modeling in cancer biology. JCO Clin Cancer Inform. (2019) 3:1–13. doi: 10.1200/CCI.18.00069

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Norton KA, Gong C, Jamalian S, Popel AS. Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes. (2019) 7:1. doi: 10.3390/pr7010037

CrossRef Full Text | Google Scholar

35. Cappuccio A, Tieri P, Castiglione F. Multiscale modelling in immunology: a review. Brief Bioinform. (2016) 17:408–18. doi: 10.1093/bib/bbv012

CrossRef Full Text | Google Scholar

36. Masoudi-Nejad A, Bidkhori G, Hosseini Ashtiani S, Najafi A, Bozorgmehr JH, Wang E, et al. Cancer systems biology and modeling: microscopic scale and multiscale approaches. Semin Cancer Biol. (2015) 30:60–9. doi: 10.1016/j.semcancer.2014.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Rejniak KA, McCawley LJ. Current trends in mathematical modeling of tumor-microenvironment interactions: a survey of tools and applications. Exp Biol Med. (2010) 235:411–23. doi: 10.1258/ebm.2009.009230

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ramanujan S, Pluen A, McKee TD, Brown EB, Boucher Y, Jain RK. Diffusion and convection in collagen gels: implications for transport in the tumor interstitium. Biophys J. (2002) 83:1650–60. doi: 10.1016/S0006-3495(02)73933-7

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Kihara T, Ito J, Miyake J. Measurement of biomolecular diffusion in extracellular matrix condensed by fibroblasts using fluorescence correlation spectroscopy. PLoS ONE. (2013) 8:e82382. doi: 10.1371/journal.pone.0082382

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Meulemans A, Paycha F, Hannoun P, Vulpillat M. Measurement and clinical and pharmacokinetic implications of diffusion coefficients of antibiotics in tissues. Antimicrob Agents Chemother. (1989) 33:1286–90. doi: 10.1128/AAC.33.8.1286

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Netti PA, Berk DA, Swartz MA, Grodzinsky AJ, Jain RK. Role of extracellular matrix assembly in interstitial transport in solid tumors. Cancer Res. (2000) 60:2497–503.

PubMed Abstract | Google Scholar

42. Sniekers YH, van Donkelaar CC. Determining diffusion coefficients in inhomogeneous tissues using fluorescence recovery after photobleaching. Biophys J. (2005) 89:1302–7. doi: 10.1529/biophysj.104.053652

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Anderson JD. Computational Fluid Dynamics: The Basics With Applications. New York, NY: McGraw-Hill, Inc. (1995).

Google Scholar

44. Anderson DA, Tannehill JC, Pletcher RH, Munipalli M, Shankar V. Computational Fluid Mechanics and Heat Transfer. 4th Edition. Boca Raton, FL: CRC Press (2021). doi: 10.1201/9781351124027

CrossRef Full Text | Google Scholar

45. Hirsch C. Numerical Computation of Internal and External Flows, 2nd Edn. Burlington, MA: Butterworth-Heinemann (2007).

Google Scholar

46. Marshall LE, Koomullil RP, Frost AR, Berry, Computational JL. Computational and experimental analysis of fluid transport through three-dimensional collagen-matrigel hydrogels. Ann Biomed Eng. (2017) 45:1027–38. doi: 10.1007/s10439-016-1748-6

PubMed Abstract | CrossRef Full Text | Google Scholar

47. ANSYS. ANSYS Fluent Theory Guide. Canonsburg, PA: ANSYS, Inc. (2020).

Google Scholar

48. Biorender. Create Professional Science Figures in Minutes. (2020). Available online at: https://Biorender.com/ (accessed September 20, 2020).

Google Scholar

49. Baptista AM, Adams EE, Gresho P. Benchmarks for the transport equation: the convection-diffusion forum and beyond, in quantitative skill assessment for coastal ocean models. Coastal Estuarine Stud. (1995) 6:241–68. doi: 10.1029/CE047p0241

CrossRef Full Text | Google Scholar

50. Okubo A, Karweit MJ. Diffusion from a continuous source in a uniform shear flow1. Limnol Oceanography. (1969) 14:514–20. doi: 10.4319/lo.1969.14.4.0514

CrossRef Full Text | Google Scholar

51. Hough KP, Wilson LS, Trevor JL, Strenkowski JG, Maina N, Kim YI, et al. Unique lipid signatures of extracellular vesicles from the airways of asthmatics. Sci Rep. (2018) 8:10340. doi: 10.1038/s41598-018-28655-9

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wang Y, Schafer CC, Hough KP, Tousif S, Duncan SR, Kearney JF, et al. Myeloid-derived suppressor cells impair B cell responses in lung cancer through IL-7 and STAT5. J Immunol. (2018) 201:278–95. doi: 10.4049/jimmunol.1701069

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Schafer CC, Wang Y, Hough KP, Sawant A, Grant SC, Thannickal VJ, et al. Indoleamine 2,3-dioxygenase regulates anti-tumor immunity in lung cancer by metabolic reprogramming of immune cells in the tumor microenvironment. Oncotarget. (2016) 7:75407–24. doi: 10.18632/oncotarget.12249

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Sawant A, Schafer CC, Jin TH, Zmijewski J, Tse HM, Roth J, et al. Enhancement of antitumor immunity in lung cancer by targeting myeloid-derived suppressor cell pathways. Cancer Res. (2013) 73:6609–20. doi: 10.1158/0008-5472.CAN-13-0987

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Sawant A, Deshane J, Jules J, Lee CM, Harris BA, Feng X, et al. Myeloid-derived suppressor cells function as novel osteoclast progenitors enhancing bone loss in breast cancer. Cancer Res. (2013) 73:672–82. doi: 10.1158/0008-5472.CAN-12-2202

PubMed Abstract | CrossRef Full Text | Google Scholar