# 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] Xuhui Peng, Rangrang Zhang. Approximations of stochastic 3D tamed Navier-Stokes equations. Communications on Pure & Applied Analysis, 2020, 19 (12) : 5337-5365. doi: 10.3934/cpaa.2020241 [2] Yi An, Bo Li, Lei Wang, Chao Zhang, Xiaoli Zhou. Calibration of a 3D laser rangefinder and a camera based on optimization solution. Journal of Industrial & Management Optimization, 2021, 17 (1) : 427-445. doi: 10.3934/jimo.2019119 [3] Balázs Kósa, Karol Mikula, Markjoe Olunna Uba, Antonia Weberling, Neophytos Christodoulou, Magdalena Zernicka-Goetz. 3D image segmentation supported by a point cloud. Discrete & Continuous Dynamical Systems - S, 2021, 14 (3) : 971-985. doi: 10.3934/dcdss.2020351 [4] Xin-Guang Yang, Lu Li, Xingjie Yan, Ling Ding. The structure and stability of pullback attractors for 3D Brinkman-Forchheimer equation with delay. Electronic Research Archive, 2020, 28 (4) : 1395-1418. doi: 10.3934/era.2020074 [5] Cung The Anh, Dang Thi Phuong Thanh, Nguyen Duong Toan. Uniform attractors of 3D Navier-Stokes-Voigt equations with memory and singularly oscillating external forces. Evolution Equations & Control Theory, 2021, 10 (1) : 1-23. doi: 10.3934/eect.2020039 [6] Xiaopeng Zhao, Yong Zhou. Well-posedness and decay of solutions to 3D generalized Navier-Stokes equations. Discrete & Continuous Dynamical Systems - B, 2021, 26 (2) : 795-813. doi: 10.3934/dcdsb.2020142 [7] Xiaoli Lu, Pengzhan Huang, Yinnian He. Fully discrete finite element approximation of the 2D/3D unsteady incompressible magnetohydrodynamic-Voigt regularization flows. Discrete & Continuous Dynamical Systems - B, 2021, 26 (2) : 815-845. doi: 10.3934/dcdsb.2020143 [8] Yang Liu. Global existence and exponential decay of strong solutions to the cauchy problem of 3D density-dependent Navier-Stokes equations with vacuum. Discrete & Continuous Dynamical Systems - B, 2021, 26 (3) : 1291-1303. doi: 10.3934/dcdsb.2020163 [9] Laurence Cherfils, Stefania Gatti, Alain Miranville, Rémy Guillevin. Analysis of a model for tumor growth and lactate exchanges in a glioma. Discrete & Continuous Dynamical Systems - S, 2020  doi: 10.3934/dcdss.2020457 [10] Duy Phan. Approximate controllability for Navier–Stokes equations in $\rm3D$ cylinders under Lions boundary conditions by an explicit saturating set. Evolution Equations & Control Theory, 2021, 10 (1) : 199-227. doi: 10.3934/eect.2020062 [11] Thomas Y. Hou, Dong Liang. Multiscale analysis for convection dominated transport equations. Discrete & Continuous Dynamical Systems - A, 2009, 23 (1&2) : 281-298. doi: 10.3934/dcds.2009.23.281 [12] Tin Phan, Bruce Pell, Amy E. Kendig, Elizabeth T. Borer, Yang Kuang. Rich dynamics of a simple delay host-pathogen model of cell-to-cell infection for plant virus. Discrete & Continuous Dynamical Systems - B, 2021, 26 (1) : 515-539. doi: 10.3934/dcdsb.2020261 [13] Olivier Pironneau, Alexei Lozinski, Alain Perronnet, Frédéric Hecht. Numerical zoom for multiscale problems with an application to flows through porous media. Discrete & Continuous Dynamical Systems - A, 2009, 23 (1&2) : 265-280. doi: 10.3934/dcds.2009.23.265 [14] Rajendra K C Khatri, Brendan J Caseria, Yifei Lou, Guanghua Xiao, Yan Cao. Automatic extraction of cell nuclei using dilated convolutional network. Inverse Problems & Imaging, 2021, 15 (1) : 27-40. doi: 10.3934/ipi.2020049 [15] Jianquan Li, Xin Xie, Dian Zhang, Jia Li, Xiaolin Lin. Qualitative analysis of a simple tumor-immune system with time delay of tumor action. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020341 [16] Yuxin Zhang. The spatially heterogeneous diffusive rabies model and its shadow system. Discrete & Continuous Dynamical Systems - B, 2020  doi: 10.3934/dcdsb.2020357 [17] Philippe Laurençot, Christoph Walker. Variational solutions to an evolution model for MEMS with heterogeneous dielectric properties. Discrete & Continuous Dynamical Systems - S, 2021, 14 (2) : 677-694. doi: 10.3934/dcdss.2020360 [18] Shigui Ruan. Nonlinear dynamics in tumor-immune system interaction models with delays. Discrete & Continuous Dynamical Systems - B, 2021, 26 (1) : 541-602. doi: 10.3934/dcdsb.2020282 [19] Shujing Shi, Jicai Huang, Yang Kuang. Global dynamics in a tumor-immune model with an immune checkpoint inhibitor. Discrete & Continuous Dynamical Systems - B, 2021, 26 (2) : 1149-1170. doi: 10.3934/dcdsb.2020157 [20] Qiang Long, Xue Wu, Changzhi Wu. Non-dominated sorting methods for multi-objective optimization: Review and numerical comparison. Journal of Industrial & Management Optimization, 2021, 17 (2) : 1001-1023. doi: 10.3934/jimo.2020009

2018 Impact Factor: 1.313