Tumor Morphology

Introduction

How to illustrate the evolution of a avascular tumor? How does it transform to a malignant tumors.

It is estimated that the new cancer cases and the deaths in 2022 are 1.9 million new cancer cases diagnosed and 609,360 cancer deaths in the United States (data from the American cancer society). The high mortality rate triggers scientists' attention. For oncologists, they develop new radiation therapy and target therapy to cure patients. However, some therapies have been proven invalid. For instance, in the previous envisage, tumors obtain nutrients by angiogenic. Thus, it is predicted that anti-angiogenic therapy provides an opportunity to eliminate cancer cells. However, it has been observed that anti-angiogenic therapy does not eliminate cancer cells but promotes tumor fragmentation and invasion [1, 2]. It seems that the growth of a tumor is much more complicated than we thought. Therefore, to better understand the surface instability of tumors, mathematical models from different perspectives, such as pressure-driven cell motion and reaction-diffusion models, are developed. In our research, we build a phenomenological model to understand how tumor growth develops invasive patterns.

To build a phenomenological model to illustrate the behavior of tumor cells, first of all, we must list all critical factors that we expect to influence tumor growth. By denoting the tumor cell density and nutrient density as $n$ and $u$, we know the death rate of tumor cells is proportional to its density, say $\gamma n$. Also, we know tumor consumes nutrient for growth. Typically, in the previous model, most of scientists assume the growth rate is linearly proportional to nutrient density and tumor cell density $n\times u$. However, it has been shown in EMT6/Ro mouse mammary tumor cells that the increase of the cell proliferation rate slows down as the nutrient density increases due to limited consumption and metabolism, and eventually, the proliferation rate saturates to a constant rate [3, 4].Therefore, we assume the growth rate to be a saturating function when nutrient density is high and use a shifted hyperbolic tangent to approximate it. That is, $$ \text{growth rate} = a\cdot n\cdot A(u), $$ where $$ A(u) = \tanh\left(\frac{u-\bar{u}}{W}\right) - A_0. $$ The coefficient $A_0$ adjusts the function since we expect a zero growth rate when there is no nutrient. The constants $\bar{u}$ and $W$ further characterize the threshold and sensitivity for tumor growth. Since the nutrient consumption rate is proportional to the tumor growth rate, the consumption rate is assumed to be $b \cdot n\cdot A(u)$, where $b$ is the proportionality coefficient.

Besides, according to previous research, the motion of tumor cells is known to be driven by a pressure gradient induced by the inner necrosis core and the outer cell-free space[1, 5]. The velocity of tumor cells obey Darcy's law $\vec{v} = -\mu\nabla P$, where $\mu$ is the mobility of tumors and $P$ is the pressure of cells. By considering the solid properties of tumor cells, the relation between pressure and density should be $dP = Kdn/n$, where $K$ is the bulk modulus of tumor cells. Thus, the velocity can be rewritten as $\vec{v} = - (\mu K/n)\nabla n$, and the flux is therefore $\nabla \cdot (n\vec{v}) = \mu K \nabla^2 n$ by assuming a homogeneous bulk modulus. The form is similar to the diffusion $D_n \nabla^2 n$. In fact, if we calculated the estimated $\mu K$, the value ranges from $1.03 \times 10^{-8}$ to $1.9 \times 10^{-8} cm^2/s$. [5, 6-9]In general, the diffusion coefficient of tumors ranges from $10^{-9$}$ to $10^{-7} cm^2/s$ depending on the type of tumors [10-12]. In our research, we apply diffusion to illustrate the motion of tumor cells. Since the motion of nutrient particles is Brownian, the motion of nutrient particles is also diffusive.

Based on the properties mentioned above, the equation of motion of tumor cells and nutrients can be written as $$ \frac{\partial n}{\partial t} = D_n \nabla^2 n + anA(u) - \gamma n - Mn(n-\alpha)(n-\beta) \\ \frac{\partial u}{\partial t} = D_u \nabla^2 u - bnA(u), $$ where the final term in the tumor evolution equation describes the difficulty of generating a tumor from a vanishing state $0$ to existing state $\beta$, the coefficient $M$ illustrates the difficulty, and the coefficient $\alpha$ illustrates the difference between the instability of the two state. For simplicity, we assume $\alpha = \beta/2$. With these equations, we can simulate the evolution of avascular tumor cells by assuming a fixed nutrient boundary condition at infinity. That is, $u(\pm\infty) = C_0$. Though it is general to simulate a three-dimensional tumor, if we focus on a symmetrical tumor, it is not necessary. By focusing on a small part on the surface of a symmetrical tumor, the system becomes a planer interface. Assuming the growth direction at $x$ axis, since the remaining axes are isotropic, we can only focus on one axis, say $y$. Moreover, since we focus on how a benign tumor with a regular shape transforms into an irregular malignant tumor, a flat surface that is homogeneous along $y$ direction is expected. Thus, we can simulate the system only on $x$ axis. The simulation is far more complicated than the traffic mode, and I remain the description on the Simulation Process page. The simulation results are shown in Fig. 1. Apparently, the tumor will grow initially until the nutrient is insufficient to maintain the tumor growth. Then, the tumor will expand to pursue more nutrients and start to die out in the nutrient-deficient region, which is generated because nutrient does not have enough time to diffuse, and outer tumors also consume a lot of nutrients. Eventually, the tumors generate two surfaces with a central necrosis core.

Fig. 1 - The evolution of tumor cells and nutrient density for $(C_0, a, \gamma, M, D_u, \bar{u}) = (8, 0.5, 0.2, 2.1, 10, 2)$

To investigate the instability of the tumor surface, assuming the flat, steady solution as $n_0$ and $u_0$, we give a small perturbation $n_1$ and $u_1$ with wave number $k$ to the flat surface. The largest eigenvalue of perturbation will determine the instability of the tumor surface. To obtain the eigenvalue and its corresponding eigenfunction, we simulate the perturbation $n_1$ and $u_1$until they are shape-preserving. The final shape-preserving eigenfunction is shown in Fig. 2.

Fig. 2 - Shape preserving eigenfunction. The coefficients used in the simulation are $(C_0, a, \gamma, M, D_u, \bar{u}) = (6, 0.5, 0.2, 2.1, 10, 2)$ for $k=0.0042$

By giving different wave numbers $k$ to the system we can obtain the dispersion relation and find the stability of the tumor surface. For positive growth rates, the tumor surface is unstable, an applied perturbation will grow, and the tumor will start to split. Fig. 3 shows the dispersion curve. Interestingly, we can find that all growth rates are negative when the boundary nutrient becomes larger.

Fig. 3 - Dispersion relation for various boundary nutrient densities ranging from $C_0 = 6$ to $C_0 = 9$. The coefficients used in the simulation are $(a, \gamma, M, D_u, \bar{u}) = (0.5, 0.2, 2.1, 10, 2)$

To understand the reason, we need to know that there are wo mechanisms that drive the instability of the system. First, the term $Mn(n-\alpha)(n-\beta)$ illustrates the boundary crossing from no tumor to existing tumor. The place occurring is on the surface, which implies that this term describes the surface energy of the system. Such a description is commonly used in many solid systems. The surface tension induced by this effect always smoothens the system, especially when the curvature is large, which can be observed when the wave number of perturbation is large. In contrast, the peaks and troughs on the perturbation will live in different nutrient environments. Peaks live in a nutritious environment, while troughs live in the nutrient-deficient region, which means the growth rates are different along peaks and troughs. This effect will enlarge the perturbations. The competition between the two effects determines the instability of the tumor surface. However, since the growth rate saturates in the nutritious region, when the nutrient is high, the effect of the growth rate difference decreases, and the system becomes surface tension dominant. Therefore, the system remains flat. The phase for different boundary nutrient values and growth coefficients aa are shown in Fig. 4.

Fig. 3 - Phase diagram for the avascular solid tumor growth. Red squares represent the phase of dying out of tumor cells. Blue circles and green triangles are breakups of tumors and steady moving tumor fronts, respectively. The coefficients used in the simulation are $(\gamma, M, D_u, \bar{u}) = (0.2, 2.1, 10, 2)$

This is an interesting result since it implies that if we plan to inhibit tumor invasion, we need to provide them with more nutrients. However, it is worth noting that tumors also grow faster when they live in a nutritious environment, and when they grow up, they might start to extend vascular. Since our research assume the tumor is avascular, further investigations are necessary to show the importance of a nutritious environment.

Further Reading

  1. V. Cristini, H. B. Frieboes, R. Gatenby, S. Caserta, M. Ferrari, and J. Sinek, Morphologic instability and cancer invasion, Clinical Cancer Research 11, 6772 (2005). (weblink)
  2. P. S. Steeg, Angiogenesis inhibitors: motivators of metastasis?, Nature Medicine 9, 822 (2003). (weblink)
  3. J. J. Casciari, S. V. Sotirchos, and R. M. Sutherland, Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular ph, Journal of Cellular Physiology 151, 386 (1992). (weblink)
  4. J. P. Freyer and R. M. Sutherland, Regulation of Growth Saturation and Development of Necrosis in EMT6/Ro Multicellular Spheroids by the Glucose and Oxygen Supply1, Cancer Research 46, 3504 (1986). (weblink)
  5. H. B. Frieboes, X. Zheng, C.-H. Sun, B. Tromberg, R. Gatenby, and V. Cristini, An integrated computational/experimental model of tumor invasion, Cancer research 66, 1597—1604 (2006). (weblink)
  6. Y. Hirakawa, T. Yoshihara, M. Kamiya, I. Mimura, D. Fujikura, T. Masuda, R. Kikuchi, I. Takahashi, Y. Urano, S. Tobita, and M. Nangaku, Quantitating intracellular oxygen tension in vivo by phosphorescence lifetime measurement, Scientific Reports 5, 17838 (2015). (weblink)
  7. M. T. Islam, S. Tang, C. Liverani, S. Saha, E. Tasciotti, and R. Righetti, Non-invasive imaging of young’s modulus and poisson’s ratio in cancers in vivo, Scientific Reports 10, 7266 (2020). (weblink)
  8. D. C. Stewart, A. Rubiano, K. Dyson, and C. S. Simmons, Mechanical characterization of human brain tumors from patients and comparison to potential surgical phantoms, PLOS ONE 12, 1 (2017). (weblink)
  9. D. Chauvet, M. Imbault, L. Capelle, C. Demene, M. Mossad, C. Karachi, A.-L. Boch, J.-L. Gennisson, and M. Tanter, In vivo measurement of brain tumor elasticity using intraoperative shear wave elastography., Ultraschall in der Medizin (Stuttgart, Germany : 1980) 37, 584 (2016). (weblink)
  10. K. R. Swanson, E. C. Alvord Jr, and J. Murray, A quantitative model for differential motility of gliomas in grey and white matter, Cell proliferation 33, 317 (2000). (weblink)
  11. S. Miyamoto, K. Atsuyama, K. Ekino, and T. Shin, Estimating the diffusion coefficients of sugars using diffusion experiments in agar-gel and computer simulations, Chemical and Pharmaceutical Bulletin 66, 632 (2018). (weblink)
  12. M. Martín-Landrove, Reaction-diffusion models for glioma tumor growth, arXiv preprint arXiv:1707.09409 (2017). (weblink)