# American Institute of Mathematical Sciences

April  2018, 15(2): 361-392. doi: 10.3934/mbe.2018016

## A multiscale model for heterogeneous tumor spheroid in vitro

 Department of Mathematical Sciences, Georgia Southern University, Statesboro, GA, 30460, USA

* Corresponding author: Zhan Chen (zchen@georgiasouthern.edu).

Received  August 30, 2016 Accepted  April 21, 2017 Published  June 2017

In this paper, a novel multiscale method is proposed for the study of heterogeneous tumor spheroid growth in vitro. The entire tumor spheroid is described by an ellipsoid-based model while nutrient and other environmental factors are treated as continua. The ellipsoid-based discrete component is capable of incorporating mechanical effects and deformability, while keeping a minimum set of free variables to describe complex shape variations. Moreover, our purely cell-based description of tumor avoids the complex mutual conversion between a cell-based model and continuum model within a tumor, such as force and mass transformation. This advantage makes it highly suitable for the study of tumor spheroids in vitro whose size are normally less than 800 $μ m$ in diameter. In addition, our numerical scheme provides two computational options depending on tumor size. For a small or medium tumor spheroid, a three-dimensional (3D) numerical model can be directly applied. For a large spheroid, we suggest the use of a 3D-adapted 2D cross section configuration, which has not yet been explored in the literature, as an alternative for the theoretical investigation to bridge the gap between the 2D and 3D models. Our model and its implementations have been validated and applied to various studies given in the paper. The simulation results fit corresponding in vitro experimental observations very well.

Citation: Zhan Chen, Yuting Zou. A multiscale model for heterogeneous tumor spheroid in vitro. Mathematical Biosciences & Engineering, 2018, 15 (2) : 361-392. doi: 10.3934/mbe.2018016
##### References:

show all references

##### References:
(A): A 3D aggregate in a hanging drop culture; (B): The representation of the Kelvin and growth elements that characterize the internal rheology of each cell, modified from previous papers [45]. Note that in a hanging drop spheroid systems in vitro, the surrounding environment exerts little resistance to growth. As such, it is reasonable to assume that no external force is imposed on in silico spheroids. Here each tumor cell inside a spheroid is modeled as a 3D deformable ellipsoid with three axes a, b, c each of which is represented by a Kelvin element. In the a-axis (similar to b-and c-axis), $u_a$ is the total change of the length, $u_a^0$ and $u_a^g$ are the changes of the length in the a-axis due to the change in the passive and growth elements respectively, $f_2$ is the nonlinear spring force from the spring in parallel, $f_a$ is the magnitude of the force applied to each end, $\mu_a$ is the viscous coefficient of the dash-pot, $k_a$ is the spring constant for the spring in the Maxwell element.
Differential adhesive forces among heterotypical cells lead to various aggregation patterns. Experimental observations of the cross section of a 3D aggregate are listed in (B1), (B2) and (B3). The corresponding numerical patterns, generated by our model after T= 7 hours for an aggregate of 1021 cells, are shown in (A1), (A2) and (A3). The same sorting pattern persists afterward. In particular, the choice of parameters $\alpha_{g,g}:\alpha_{r,r}:\alpha_{g,r}=0.4:1:0.7$ (DA 2) leads to cell sorting in (A2). A cross section of the 3D configuration is shown under its 3D counterpart. Cells do not sort at all when $\alpha_{g,g}:\alpha_{r,r}:\alpha_{g,r}=1:1:1$ in (A1). With $\alpha_{g,g}:\alpha_{r,r}:\alpha_{g,r}=1:1:0.2$ (DA 1), cells separate in (A3). (B1), (B2) and (B3) are the experimental observations from Duguay et al [16] where two L-cell lines express N-cad at different levels. Line N5A expresses about 50% more N-cad than what line N2 does. An aggregate in figure (B1) does not sort, in which both the red-and green-colored cells are from line N5A after 1 day of culture. Yet a similar aggregate in figure (B2) containing a mixture of N5A (red) and N2 (green) cells segregate from one another during 1 day of culture, where higher-expressing N5A cells were completely enveloped by lower-expressing N2 cells. In figure (B3), Aggregates containing equal numbers of L cells lead to mounds of R-cad-expressing cells (red) partially capping a B-cad-expressing mass (green) after being cultured in suspension for 2 days.
(A) 3D Simulation results of the same engulfment pattern by fragment fusion and sorting shown in a cross section of a 3D aggregate; (B) in vitro observations. In both cases, we set $\alpha_{g,g}:\alpha_{r,r}:\alpha_{g,r}=0.4:1:0.7$. In an aggregate starting with intermixed cells, cells sort by the coalescence of smaller islands to form larger ones (sorting); If two tissues have initial contact, green cells gradually spread over red ones and eventually envelop them (fragment fusion). Our in silico results reproduced the in vitro observations (B) which were taken from Foty's review [24,75] where zebrafish ectoderm and mescendoderm tissues were mixed together or contacted each other. The system reached a stable configuration after 16 h, as the ectoderm occupied the internal position.
Compression forces on cell sorting. In all these simulations, we set $\alpha_{g,g}:\alpha_{r,r}:\alpha_{g,r}=0.4:1:0.7$. Ratio S is calculated by the ratio of the total number of lighter adhesive cells over the total number of the cells in the outer part of the smallest rectangle solid containing the aggregate
Growth of the tumor spheroid with a pre-existing necrotic core based on the 2D cross section configuration. (a) the oxygen profile which is described in percentage; (b) initial configuration; (c) intermediate state (T= 24 hours); (d): final configuration (T= 48 hours) where green (or blue) is for proliferating cells, red is for quiescent cells and black is for the necrotic core. The unit of color bar in the oxygen is in percentage. One percentage is equal to 0.013 mM. The spatial unit is per 10 $\mu m$.
Plots of the stabilized viable rim (defined by one half of the difference between the diameter of a tumor and the diameter of its necrotic core) in two cases: (a) the evolution of the viable rim inside the tumor spheroid with a pre-existing necrotic core. (b) the evolution of the variable rim inside the tumor spheroid which consists only of initial proliferating cells. Radius is chosen on right plot to enable readers to observe clearly when the necrotic core emerges and how it evolves. Data point at t=0 is not measured since it takes time for an initialized spheroid system to approach a mechanical quasi-equilibrium. In addition, after the stabilization from t=2000 minutes, the same pattern as (a) was observed. But there are slight fluctuations in the simulation due to cellular random walks.
Growth of the tumor spheroid without a pre-existing necrotic core based on the 2D cross section configuration. (a) initial configuration, (b) configuration after T= 44 hours. The spatial unit is per 10 $\mu m$ here.
(A) Evolution of the frequency of labelled cells in a growing spheroid at three different elapsed time points (2 hours, 24 hours and 48 hours); (B) Evolution of the frequency of labelled cells in a non-growing spheroid. Homotypic labeled cells are initially adhered to the surface of the spheroids. The number of those cells are recorded in different depths of the spheroids.
First row for a non-growing spheroid where one type of cells experience chemotaxis: (a) initial configuration, (b) final distribution (T= 48 hours). Second row for a growing spheroid where one type of cells experience chemotaxis: (c) initial configuration, (d) final distribution (T= 48 hours). Active forces for both random motion and chemotaxis are 8 nN. The spatial unit is per 10 $\mu m$ here.
The pathway of cell sorting in a 2D cross section where $\alpha_{r,r}:\alpha_{g,g}:\alpha_{g,r}=0.4:1.5:0.7$. It can be seen that small islands of green cells fuse into large ones.
An illustration of passive force using Equation 12 Here two standard spherical cells of 10 $\mu m$ diameter are used to carry out the calculations
Parameters for the cell-based component of the model.
 Parameter Description Value Dimensionless in coding Refs. Adhesion parameters $\mu_{cell}$ cell-cell adhesiveness 27.0 dyn s/cm 450 [11,45] $\mu_s$ cell-substrate 27.0 dyn s/cm 450 [11,45] adhesiveness $\mu_{f}$ fluid viscosity 2.7 dyn s/cm 450 [11,45] Rheological parameters $c^+$ growth function 5.16089$\times 10^{-9}$ mm/(min. nN) 5.16089$\times 10^{-9}$ [45] $\sigma^+$ growth function 800 nN 800 [45] $\sigma^-$ growth function -4 nN -4 [45] $\alpha$ growth function 0.0 nN 0.0 [45] $k_a$ standard solid 163.8 dyn/cm 163800 [11,45,84] $k_2$ standard solid 147.5 dyn/cm, 147500 [11,45,84] $\mu_a$ standard solid 123 dyn min/cm 123000 [11,45,84] $f_a$ active force 10 nN 10 in this work
 Parameter Description Value Dimensionless in coding Refs. Adhesion parameters $\mu_{cell}$ cell-cell adhesiveness 27.0 dyn s/cm 450 [11,45] $\mu_s$ cell-substrate 27.0 dyn s/cm 450 [11,45] adhesiveness $\mu_{f}$ fluid viscosity 2.7 dyn s/cm 450 [11,45] Rheological parameters $c^+$ growth function 5.16089$\times 10^{-9}$ mm/(min. nN) 5.16089$\times 10^{-9}$ [45] $\sigma^+$ growth function 800 nN 800 [45] $\sigma^-$ growth function -4 nN -4 [45] $\alpha$ growth function 0.0 nN 0.0 [45] $k_a$ standard solid 163.8 dyn/cm 163800 [11,45,84] $k_2$ standard solid 147.5 dyn/cm, 147500 [11,45,84] $\mu_a$ standard solid 123 dyn min/cm 123000 [11,45,84] $f_a$ active force 10 nN 10 in this work
Sorting in the presence of differential adhesion. Various ratios of cohesion can lead to cell sorting in two types of cells, distinguished by green and red in this work. Here $\alpha_{g,g}$, $\alpha_{r,r}$, $\alpha_{g,r}$ represent relative adhesive strengths between like and unlike cells (green and green, red and red, or green and red respectively)
 Index $\alpha_{g,g}$ $\alpha_{r,r}$ $\alpha_{g,r}$ Sorting results 1 0.4 1 0.7 green envelops red 2 0.6 1 0.8 green envelops red 3 0.8 1 0.9 green envelops red 4 1 1 1 Not sorting 5 1 1 0.2 green and red separate
 Index $\alpha_{g,g}$ $\alpha_{r,r}$ $\alpha_{g,r}$ Sorting results 1 0.4 1 0.7 green envelops red 2 0.6 1 0.8 green envelops red 3 0.8 1 0.9 green envelops red 4 1 1 1 Not sorting 5 1 1 0.2 green and red separate
Parameters used in the reaction-diffusion component of the model. We use the cell average packing density carried out $2.01\times10^8\;\;\mbox{cells}/cm^3$ in Casciari $et al$ . [8] to convert uptake parameters $A_{O_2},A_{gl},B_{O_2},B_{gl}$ in this table to rates per unit volume.
 P Description Value Dimensionless in coding Refs. Diffusion Coefficients of oxygen in each region $D_o^c$ cell based region $1.82\times 10^{-5}\;cm^2/s$ 6.552 [45] $D_o^q$ continuum region $2.15\times 10^{-6}\; cm^2/s$ 7.74 Diffusion Coefficients of glucose in each region $D_g^p$ cell based region $3.0\times 10^{-6}\; cm^2/s$ 1.08 this work $D_g^q$ continuum region $6.46\times 10^{-6}\; cm^2/s$ 2.3256 [45] Coefficients in Uptake Functions $A_{O_2}$ oxygen uptake $1.0642\times 10^{-16}\; \frac{mol}{cell\cdot s}$ 2.01014 [9,45] $B_{O_2}$ oxygen uptake $6.0202\times 10^{-17}\; \frac{mol\cdot mM}{cell\cdot s}$ 0.0497 [8,9,45] $A_{gl}$ glucose uptake $1.0642\times 10^{-16}\; \frac{mol}{cell\cdot s}$ 2.01014 [8,9,45] $B_{gl}$ glucose uptake $1.7879\times 10^{-17}\; \frac{mol\cdot mM}{cell\cdot s}$ 0.0107 [8,25,45] $k_{O_2}$ critical oxygen concentration $4.640\times 10^{-3}\; mM$ $1.856\times 10^{-4}$ [8,45] $k_{gl}$ critical glucose concentration $4.0\times 10^{-2}\; mM$ $1.6\times 10^{-3}$ [8,45] $n_{O_2}$ oxygen uptake $0.55\; mM$ $2.2 \times 10 ^{-2}$ [25,45] $n_{gl}$ glucose uptake $0.04\; mM$ $1.6 \times 10 ^{-3}$ [25,45]
 P Description Value Dimensionless in coding Refs. Diffusion Coefficients of oxygen in each region $D_o^c$ cell based region $1.82\times 10^{-5}\;cm^2/s$ 6.552 [45] $D_o^q$ continuum region $2.15\times 10^{-6}\; cm^2/s$ 7.74 Diffusion Coefficients of glucose in each region $D_g^p$ cell based region $3.0\times 10^{-6}\; cm^2/s$ 1.08 this work $D_g^q$ continuum region $6.46\times 10^{-6}\; cm^2/s$ 2.3256 [45] Coefficients in Uptake Functions $A_{O_2}$ oxygen uptake $1.0642\times 10^{-16}\; \frac{mol}{cell\cdot s}$ 2.01014 [9,45] $B_{O_2}$ oxygen uptake $6.0202\times 10^{-17}\; \frac{mol\cdot mM}{cell\cdot s}$ 0.0497 [8,9,45] $A_{gl}$ glucose uptake $1.0642\times 10^{-16}\; \frac{mol}{cell\cdot s}$ 2.01014 [8,9,45] $B_{gl}$ glucose uptake $1.7879\times 10^{-17}\; \frac{mol\cdot mM}{cell\cdot s}$ 0.0107 [8,25,45] $k_{O_2}$ critical oxygen concentration $4.640\times 10^{-3}\; mM$ $1.856\times 10^{-4}$ [8,45] $k_{gl}$ critical glucose concentration $4.0\times 10^{-2}\; mM$ $1.6\times 10^{-3}$ [8,45] $n_{O_2}$ oxygen uptake $0.55\; mM$ $2.2 \times 10 ^{-2}$ [25,45] $n_{gl}$ glucose uptake $0.04\; mM$ $1.6 \times 10 ^{-3}$ [25,45]
 [1] Gülnihal Meral, Christian Stinner, Christina Surulescu. On a multiscale model involving cell contractivity and its effects on tumor invasion. Discrete & Continuous Dynamical Systems - B, 2015, 20 (1) : 189-213. doi: 10.3934/dcdsb.2015.20.189 [2] Yangjin Kim, Hans G. Othmer. Hybrid models of cell and tissue dynamics in tumor growth. Mathematical Biosciences & Engineering, 2015, 12 (6) : 1141-1156. doi: 10.3934/mbe.2015.12.1141 [3] Thomas Y. Hou, Zuoqiang Shi. Dynamic growth estimates of maximum vorticity for 3D incompressible Euler equations and the SQG model. Discrete & Continuous Dynamical Systems - A, 2012, 32 (5) : 1449-1463. doi: 10.3934/dcds.2012.32.1449 [4] Ahuod Alsheri, Ebraheem O. Alzahrani, Asim Asiri, Mohamed M. El-Dessoky, Yang Kuang. Tumor growth dynamics with nutrient limitation and cell proliferation time delay. Discrete & Continuous Dynamical Systems - B, 2017, 22 (10) : 3771-3782. doi: 10.3934/dcdsb.2017189 [5] Yong Zhou. Remarks on regularities for the 3D MHD equations. Discrete & Continuous Dynamical Systems - A, 2005, 12 (5) : 881-886. doi: 10.3934/dcds.2005.12.881 [6] Hyeong-Ohk Bae, Bum Ja Jin. Estimates of the wake for the 3D Oseen equations. Discrete & Continuous Dynamical Systems - B, 2008, 10 (1) : 1-18. doi: 10.3934/dcdsb.2008.10.1 [7] Indranil SenGupta, Weisheng Jiang, Bo Sun, Maria Christina Mariani. Superradiance problem in a 3D annular domain. Conference Publications, 2011, 2011 (Special) : 1309-1318. doi: 10.3934/proc.2011.2011.1309 [8] Giovanny Guerrero, José Antonio Langa, Antonio Suárez. Biodiversity and vulnerability in a 3D mutualistic system. Discrete & Continuous Dynamical Systems - A, 2014, 34 (10) : 4107-4126. doi: 10.3934/dcds.2014.34.4107 [9] Sadek Gala. A new regularity criterion for the 3D MHD equations in $R^3$. Communications on Pure & Applied Analysis, 2012, 11 (3) : 973-980. doi: 10.3934/cpaa.2012.11.973 [10] Jiahong Wu. Regularity results for weak solutions of the 3D MHD equations. Discrete & Continuous Dynamical Systems - A, 2004, 10 (1&2) : 543-556. doi: 10.3934/dcds.2004.10.543 [11] Gabriel Deugoue. Approximation of the trajectory attractor of the 3D MHD System. Communications on Pure & Applied Analysis, 2013, 12 (5) : 2119-2144. doi: 10.3934/cpaa.2013.12.2119 [12] Alp Eden, Varga K. Kalantarov. 3D convective Cahn--Hilliard equation. Communications on Pure & Applied Analysis, 2007, 6 (4) : 1075-1086. doi: 10.3934/cpaa.2007.6.1075 [13] Chongsheng Cao. Sufficient conditions for the regularity to the 3D Navier-Stokes equations. Discrete & Continuous Dynamical Systems - A, 2010, 26 (4) : 1141-1151. doi: 10.3934/dcds.2010.26.1141 [14] Ning Ju. The global attractor for the solutions to the 3D viscous primitive equations. Discrete & Continuous Dynamical Systems - A, 2007, 17 (1) : 159-179. doi: 10.3934/dcds.2007.17.159 [15] Tomás Caraballo, Antonio M. Márquez-Durán, José Real. Pullback and forward attractors for a 3D LANS$-\alpha$ model with delay. Discrete & Continuous Dynamical Systems - A, 2006, 15 (2) : 559-578. doi: 10.3934/dcds.2006.15.559 [16] Jianqing Chen. Best constant of 3D Anisotropic Sobolev inequality and its applications. Communications on Pure & Applied Analysis, 2010, 9 (3) : 655-666. doi: 10.3934/cpaa.2010.9.655 [17] Ming Lu, Yi Du, Zheng-An Yao. Blow-up phenomena for the 3D compressible MHD equations. Discrete & Continuous Dynamical Systems - A, 2012, 32 (5) : 1835-1855. doi: 10.3934/dcds.2012.32.1835 [18] Felipe Linares, Jean-Claude Saut. The Cauchy problem for the 3D Zakharov-Kuznetsov equation. Discrete & Continuous Dynamical Systems - A, 2009, 24 (2) : 547-565. doi: 10.3934/dcds.2009.24.547 [19] Rafel Prohens, Antonio E. Teruel. Canard trajectories in 3D piecewise linear systems. Discrete & Continuous Dynamical Systems - A, 2013, 33 (10) : 4595-4611. doi: 10.3934/dcds.2013.33.4595 [20] Makram Hamouda, Chang-Yeol Jung, Roger Temam. Asymptotic analysis for the 3D primitive equations in a channel. Discrete & Continuous Dynamical Systems - S, 2013, 6 (2) : 401-422. doi: 10.3934/dcdss.2013.6.401

2018 Impact Factor: 1.313