Next Article in Journal
Cancer and Radiosensitivity Syndromes: Is Impaired Nuclear ATM Kinase Activity the Primum Movens?
Next Article in Special Issue
Endogenous Extracellular Matrix Regulates the Response of Osteosarcoma 3D Spheroids to Doxorubicin
Previous Article in Journal
Therapeutic Targeting of MERTK and BCL-2 in T-Cell and Early T-Precursor Acute Lymphoblastic Leukemia
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Bio-Mechanical Model of Osteosarcoma Tumor Microenvironment: A Porous Media Approach

Department of Mathematics and Statistics, University of Massachusetts Amherst, Amherst, MA 01003, USA
*
Author to whom correspondence should be addressed.
Cancers 2022, 14(24), 6143; https://doi.org/10.3390/cancers14246143
Submission received: 10 November 2022 / Revised: 5 December 2022 / Accepted: 8 December 2022 / Published: 13 December 2022
(This article belongs to the Special Issue Osteosarcoma Microenvironment)

Abstract

:

Simple Summary

Osteosarcoma is the most common type of bone cancer seen in children to young adults with poor prognosis. To find effective treatments, it is crucial to understand the mechanism of the initiation and progression of the osteosarcoma tumors. In this paper, we introduce a PDE model for the progression of osteosarcoma tumors by considering the location of different cell types, including immune and cancer cells, in the tumors. We perform several simulations using the developed model to investigate the importance and role of the different immune cells’ location in the growth of the tumors. The results show that the co-localization of macrophages and cancer cells promotes tumors’ growth.

Abstract

Osteosarcoma is the most common malignant bone tumor in children and adolescents with a poor prognosis. To describe the progression of osteosarcoma, we expanded a system of data-driven ODE from a previous study into a system of Reaction-Diffusion-Advection (RDA) equations and coupled it with Biot equations of poroelasticity to form a bio-mechanical model. The RDA system includes the spatio-temporal information of the key components of the tumor microenvironment. The Biot equations are comprised of an equation for the solid phase, which governs the movement of the solid tumor, and an equation for the fluid phase, which relates to the motion of cells. The model predicts the total number of cells and cytokines of the tumor microenvironment and simulates the tumor’s size growth. We simulated different scenarios using this model to investigate the impact of several biomedical settings on tumors’ growth. The results indicate the importance of macrophages in tumors’ growth. Particularly, we have observed a high co-localization of macrophages and cancer cells, and the concentration of tumor cells increases as the number of macrophages increases.

1. Introduction

Osteosarcoma (also called osteogenic sarcoma) is the most common type of cancer that starts in the bones. It happens most often in children, adolescents, and young adults. Approximately 800 new cases of osteosarcoma are reported each year in the U.S.; of these cases, about 400 are in children and teens [1]. Surgery, chemotherapy, radiation therapy, and targeted therapy are the types of standard treatment for osteosarcoma [2].
There are many mathematical models to understand and investigate biomedical problems such as cancer. Mainly, bio-mechanical models are developed to investigate the spatial interaction among cells and their movements in tumors using appropriate physical laws [3,4]. In one study, the authors have developed a partial differential equation (PDE) model to study the breast tumors’ progression in mice as a fluid structure because the breast tumors are mainly confined in the mammary gland [5]. However, when the tissues are porous, modeling with techniques of porous media can be more realistic [6,7,8,9,10,11]. Although the Biot equations are the critical part of poroelasticity theory, few studies have incorporated the equations for porous media (the Biot equations) into their model.
There is a limited number of mathematical models for osteosarcoma, and to the best of our knowledge, there is no PDE model for osteosarcoma. To find optimal treatments for osteosarcoma patients based on their immune profiles, an ODE model has been recently developed to study the key players in the tumors’ microenvironment and their interaction network, which are shown in Figure 1. In this model, osteosarcoma tumors have been grouped based on their immune profile. For each group, the models’ parameters that could not be found in current literature have been estimated using the tumors’ gene expression profiles [12,13].
Here, we propose a bio-mechanical multiphase model that describes the dynamics of the tumor microenvironment of osteosarcoma. The model PDE consists of two parts: the biological part and the mechanical part. The biological part is an extension of the ODE system provided in [12] to a system of RDA equations. The coupling terms from the mechanical part are the fluid velocity involved in the convection terms. The mechanical part is derived from the first principles, following the same argument that leads to Biot equations. Using the mechanical part, we model the motion of the fluid, which carries various kinds of cells and chemokines, and is an essential part of the continuity equation as fluid flow through porous media. The coupling term from the biological part is the additional strain of the solid tumor caused by the increase or reduction of cells. We investigate the reference case in the absence of influx and cases with different settings to examine the importance of considering cell locations and their movements in the mathematical modeling of the tumors’ growth.

2. Materials and Methods

2.1. Tumor Microenvironment and Interaction Network

The key players of the tumor microenvironment in osteosarcoma and their interaction network are shown in Figure 1, which lay the foundation of the ODE paper [12].
The vector, [ X ] = ( [ X 1 ] , [ X 2 ] , , [ X 14 ] ) represents the 14 key players we will investigate. The values of [ X i ] ’s are re-scaled by the corresponding steady-state abundance of cells or cytokines. The correspondence of [ X i ] ’s to the cells or molecules and the steady-state abundance values are shown in Table 1.

2.2. Biological Part of the Model

We adopt the Reaction-Diffusion-Advection equations by adding the diffusion term D c e l l Δ [ X i ] and the advection term b · ( [ X i ] κ p ) to the ODEs of each cell type. Here, D c e l l is the diffusion coefficient of cells, i = 2 , 4 , 5 , 6 , 7 , 8 , 9 , b is the advection coefficient which is set to be 1 for the sake of simplicity, and κ is the hydraulic conductivity of the solid tumor.The derivation of the form of the advection term will be explained in detail in the mechanical modeling part. We assume that necrotic cells are immotile, so we do not consider any diffusion or advection for them [14]. As for the naive T-cells and naive macrophages, since they activate mainly outside of the tumor microenvironment [15,16,17], we are only interested in their level, not their spatial distributions. We, therefore model their dynamics via ODEs.
We add the diffusion and advection terms to every ODE of cytokines. We denote the diffusion coefficient of HMGB1 (H) by D H and the diffusion coefficients of other cytokines ( I γ , μ 1 , and μ 2 ) by D c y t o .
In summary, the biological part of our model can be written as
[ X i ] t D i Δ [ X i ] + b i · ( [ X i ] κ p ) = f i , i = 1 , , 14 ,
where D i is 0 , D c e l l , D c y t o , or D H , depending on the specific equation, b i = 0 or 1 , and f i ’s model the biochemical interactions which are given in [12] (please see Appendix A). Consequently, the equations can be written in the following vector form
[ X ] t D Δ [ X ] + B · ( [ X ] κ p ) = f ,
where D = diag { D i } and B = diag { b i } are diagonal coefficient matrices.

2.3. Mechanical Part of the Model

In the chapter “Cancer Models and Their Mathematical Analysis” of the book [14], Friedman utilizes the continuity equation that generally takes the form
ρ t + · ( ρ v ) = f ,
where ρ is the density of the matter of interest, v is usually the fluid velocity, and f is the external source. In the context of cancer modeling, v represents the continuous motion of cells within the tumor. Later, Friedman specifies the meaning of v to be the velocity of the fluids that flow through porous media and applies Darcy’s law, treating the tumor tissue as a porous medium and assuming the moving cells are carried with the flow.
Here, we will go one step further and take advantage of the full poroelasticity equations that cover Darcy’s law while maintaining a variable that describes the movement of the solid body—the tumor. Therefore, we put the description of fluid flow into the theory of porous media and obtain a source for the solid/mesh movement.

2.3.1. Governing Equations

We start with presenting the stress( σ i j )-strain( e i j ) relation for linear, isotropic poroelastic material. In the spirit of Hooke’s law, the linear relation is
σ i j = K u 2 G 3 δ i j e + 2 G e i j α δ i j ζ .
Here, e i j = 1 2 u i x j + u j x i is the stress, e = e i i = u 1 x 1 + u 2 x 2 = · u is the total volumetric strain, where u is the solid displacement of the material, δ i j is the Kronecker delta, ζ is the fluid content, K u is the undrained bulk modulus, G is the shear modulus, and α is the Biot effective stress coefficient. The total stress comprises two parts, the contribution from the displacement of the solid framework and the contribution from the fluid flow. In our case, when building a multiphase model, we consider the contribution from the variation of the number of cells that constitute the solid frame of the tumor, namely, the cancer cells C, the necrotic cells N, the naive dendritic cells D n , and the dendritic cells D, since most tumors are infiltrated by dendritic cells [18].
The expression for the stress can be formally written as
stress = solid contribution + fluid contribution + biological contribution .
For the biological contribution, we apply Hooke’s law again. The spring constant is the drained bulk modulus K, and the distance is induced by the volumetric change resulting from the variations of cells’ number. The number of cells changes with respect to time in a pattern that we have already given:
[ X i ] t D i Δ [ X i ] + b i · ( [ X i ] v f ) = f i . i = 1 , , 14 .
The volume of cells that are a part of the solid tumor V is calculated through
V = i = 7 10 A X i [ X i ]
where A X i is the size of the cell X i . The variable X i when i = 7 , 8 , 9 , and 10 corresponds to naive dendritic cells, dendritic cells, cancer cells, and necrotic cells. To summarize, if we regard the tumor as isotropic, the stress-strain relation for the tumor is,
σ i j = K u 2 G 3 δ i j e + 2 G e i j α δ i j ζ K V δ i j .
This expression is verified in the paper by Roose et al. [19], where they have used the same form. Here, we have the density of cells at any time from solving the biological equations, so we can calculate the volumetric change without resorting to the other given formulae of tumor growth. We complete the set of governing equations with the following ones.
Definition of strain:
e i j = 1 2 u i x j + u j x i ,
where we define u as the solid displacement of the solid tumor, a macroscopic measure.
Constitutive equation of the fluid pressure p:
p = M ( ζ α e ) ,
where M is the Biot modulus.
Equilibrium equation (Newton’s law of motion):
σ i j x j = 0 .
Darcy’s law (flow of a fluid with constant viscosity across a solid is a function of its pressure difference and the solid properties):
q i = κ p x i .
Here, q is the fluid-specific discharge vector, defined as the volume of fluid passing through a unit area of the porous medium, per unit time, in the direction normal to the area, which has the dimension of velocity, [L/T], and κ is the hydraulic conductivity.
Continuity equation (conservation of mass):
ζ t + q i x i = 0 .
Solving the governing Equations (7)–(12) simultaneously, we obtain the mechanical part of our model: In Ω ( t )
K 2 G 3 ( · u ) G Δ u + α p + K V = 0 ,
t c 0 p + α · u κ Δ p = 0 ,
which is to be aligned with the biological part (1)
[ X i ] t D i Δ [ X i ] + b i · ( [ X i ] κ p ) = f i , i = 1 , , 14 .
The parameter values of the local dynamics f of the biological part are shown in [12] and other parameter values are shown in Table 2.

2.3.2. Initial and Boundary Conditions

For the mechanical part, for simplicity, we set the initial displacement field to zero and the interstitial fluid pressure equal to the blood pressure everywhere. This yields
u ( x , 0 ) = ( 0 , 0 ) , p ( x , 0 ) = p 0 ,
where p 0 is the average blood pressure of humans. Since we set u at time 0 as the reference, we shall reset the biological coupling term in Equation (13) as K V * = K ( V V 0 ) where V 0 is V ( t = 0 ) so that at every time, V is to be compared with V 0 . At t = 0 , it can be seen that u = 0 ; hence compatibility is ensured. The initial conditions of the biological part can be found in [12]; we have also listed them here in Table A1.
We use the no-flux boundary condition for u , in which case the tumor can grow freely while ignoring mechanical influences from surrounding tissues. We also treat the boundary of the tumor as the boundary of the normal tissues and set the pressure to be fixed as the blood pressure.
u ( x , t ) n = 0 , on Ω ,
p ( x , t ) = p 0 , on Ω ,
where n is the outward unit normal vector. To exclude translation or rotation of the domain, two constraints are added [24]
Ω u d x = 0 , Ω u × d d x = 0 ,
where d is the deformation vector.
We will apply both the no-flux boundary conditions and the Robin Boundary conditions to the biological equations. A no-flux boundary condition resembles the case when the tumors develop without interference, while the Robin boundary condition imitates the case when immune cells infiltrate. They can be written as
[ X i ] n = 0 , on Ω , or ,
[ X i ] n + α i ( [ X i ] [ X i ^ ] ) = 0 , on Ω ,
where α i is the influx rate and is only non-zero for immune cells. The quantity [ X i * ] pertains to the maximum levels of immune cells in lymph nodes and blood.

2.4. Weak Formulation and Discretization

2.4.1. Weak Formulation

We first rewrite the biological equations in the Eulerian system. Noticing that the definition of the material derivative is given by
D D t = t + v f · ,
we have
[ X i ] t D i Δ [ X i ] + b i · ( [ X i ] v f ) = [ X i ] t D i Δ [ X i ] + b i [ X i ] · v f + b i [ X i ] · v = D [ X i ] D t D i Δ [ X i ] + b i [ X i ] · v f = f i ,
where v f is the fluid velocity. Let I = ( 0 , T ] , V = ( L 2 ( I ; H 1 ( Ω ) ) ) 3 , Q = L 2 ( I ; H 1 ( Ω ) ) , X = ( L 2 ( I ; H 1 ( Ω ) ) ) 14 , and a u , v = K 2 G 3 · u , · v + G u , v , the weak formulation of the problem is: find ( u , p , [ X ] ) V × Q × X that satisfies (16) and Lagrange multipliers λ i such that for a.e. t ( 0 , T ]
a u , v + α p , v + K V * , v i = 1 3 λ i v , z i i = 1 3 ω i u , z i = 0 ,
c 0 p t , q + α ( · u ) t , q + κ p , q = 0 ,
and
D [ X ] D t , Ξ + D [ X ] , Ξ + b [ X ] κ Δ p , Ξ = ( f , Ξ ) ,
for all appropriate choices of test functions v , q and Ξ . The last two terms at the left-hand side of Equation (24) are added to exclude rigid body movements, with the following constraints of Lagrange multipliers λ i and test functions ω i corresponding to the bases z i of the space of rigid body movements:
i = 1 3 λ i v h , z i + i = 1 3 ω i u h , z i = 0 ,
where z i { ( 1 , 0 ) , ( 0 , 1 ) , ( x 2 , x 1 ) } .

2.4.2. Discretization

We only list the result of discretization in this section. The details about discretization are given in Appendix B.
Denote by V h , Q h , X h the appropriate finite dimensional subspaces of V , Q , X , respectively, and denote the time step size by k. Let u n = u ( t n ) and p n = p ( t n ) and let U n and P n be their approximations, respectively. We use Lagrange elements for the two mechanical equations. The fully discrete approximation for the mechanical problem is: find U n V h and P n Q h for n = 1 , 2 , , N , and Lagrange multipliers λ i such that
a U n , v h + P n , v h + K V * n , v h i = 1 3 λ i v h , z i i = 1 3 ω i u h , z i = 0 ,
k 1 c 0 ( P n P n 1 ) , q h + k 1 α · ( U n U n 1 ) , q h + κ P n , q h = 0 ,
for all basis functions v h and q h of the spaces V h and Q h .
After solving the mechanical problem (28) and (29), we obtain the mesh displacement vector U n . We can apply it to the domain if we want to observe the shape and size of the domain with respect to the reference frame, at any time.
We then solve the biological problem. We use a mixed finite element space enriched with bubble elements for the biological equations. The fully discrete approximation for the biological problem is: find X n Q h for n = 1 , 2 , , N , such that
k 1 [ X ] n [ X ] n 1 , Ξ h + D [ X ] n , Ξ h + b [ X ] n κ Δ p , Ξ h = ( f , Ξ h )
for all basis functions Ξ h of the space X h .

3. Results

There are three distinct groups of osteosarcoma tumors based on their immune profile [25]. The recent ODE model [12] investigates the dynamics of key players given in Table 1 for each of these three groups by finding their parameter values using tumors’ gene expression profiles. In this paper, we have used the parameter values for the biochemical interactions ( f i ’s) obtained for the first group of osteosarcoma tumors using the TARGET data with 88 samples (downloaded from the UCSC Xena web portal [26]). We simulated several scenarios to evaluate the model’s performance and gather insights about the role of the spatial interactions among key players and their movements during the tumors’ growth. We performed four main categories of simulations: the reference case Section 3.1, with a different initial profile Section 3.2, with immune cells infiltration from the boundary Section 3.3, and with external sources Section 3.4, Section 3.5, Section 3.6 and Section 3.7, as shown below.

3.1. The Reference Case

The reference case has the simplest settings. In this case, the initial conditions of the biological Equations (2) are uniform throughout the domain with no-flux boundary condition, meaning that there is no infiltration of cells. Finally, there is no alteration of biological coefficients, i.e., no treatments are modeled. Mathematically, we solve the Equations (13)–(15) with the boundary conditions:
t ( x , t ) · n = 0 , on Ω ,
p ( x , t ) = p 0 , on Ω ,
[ X ] n = 0 , on Ω ,
and initial conditions:
u ( x , 0 ) = ( 0 , 0 ) ,
p ( x , 0 ) = p 0 ,
[ X ] ( x , 0 ) = C i ,
where C i , for i = 1 , 2 , or 3 is the vector of initial condition for 3 clusters, shown in Table A1. The domain Ω is initially a circle with a diameter 0.01 (m).
As expected, the number of cells is the same everywhere throughout the domain at any time for all variables because the initial values of variables are uniform throughout the entire domain. The values of C and D in the entire domain are shown in Figure 2a,d. The perfect agreement between the results of the reference case (Figure 2a,e and the ones obtained from the ODE model (green curves in Figure 2c,f) provide validation for the PDE model and its implementation.
The size change of the domain is partly displayed in Figure 3. At t = 10 , the number of cells V decreases, implying that the size should decrease accordingly. The inward pointing displacement vector field u shown in Figure 3b suggests the same trend. An opposite example is given in Figure 3c, for t = 400 . The domain at t = 10 , 200 and 400 are shown in the same frame in Figure 3d–f, respectively.
The solution of the reference case is essential because it can be used to validate the model and compare the results of different scenarios with the reference case. We calculate the total number of cells that are a part of the solid tumor: cancer cells, necrotic cells, dendritic and naive dendritic cells. We aim to simulate the tumor’s size changes, as indicated by the number of cells, by applying the solid displacement vector u to the domain. Because the domain is selected to be the tumor or part of it, the solid displacement vector u tells us how each point on the solid tumor moves with respect to the reference frame at any time.
We use the initial total number of cells V 0 and the initial diameter d 0 = 0.01 ( m ) as references. If, at one time, the total number of cells V reaches n V 0 , we have an intuitive but crude estimate of the tumor size at that time, A = n A 0 = π 4 × 10 4 n ( m 2 ) , and an estimate of the diameter of the tumor d = n d 0 (since we are working on 2D cases). For example, by the time the total number of cells becomes 4 V 0 , we expect the diameter to reach approximately 0.02   ( m ) . In fact, the size can be bigger because of the remnants of dead cells inside the tumor or stronger angiogenesis, or it can be smaller if the squeezing effect dominates. In [27], the authors have built a generalized linear model that connects the tumor cell number with cell diameter and tumor diameter after studying 38 tumor samples with R 2 = 0.92 . The relation is
No . of cells / colony = 2.40 ( colony diameter ) 2.378 ( cell diameter ) 2.804 .
Using this formula, we obtain another estimate of the tumor’s diameter (Table 3). Figure 4 provides a visualized comparison between different estimates. The number of cells decreases in the first 17 days due to a sharp reduction of necrotic cells, then increase rapidly, and finally remains unchanged at about 1100 days. We, therefore, can estimate the evolution of a biological feature, the size of the tumor, using a calculated mechanical quantity, the displacement vector.

3.2. The Case of a High Level of M and T c in the Middle of the Tumor at the Initial Time

This is to imitate the condition of a high initial biological activity in the middle of the tumor. The initial concentration of M and T c in the center is set to be 11 times higher than the boundary, as shown in Figure 5a,c. The distribution of T c over the domain rapidly changes, as we can see from Figure 5d at t = 2.5 . The concentration of T c becomes higher on the boundary, and the overall concentration greatly decreases. However, M maintains the “higher concentration in the middle” profile throughout the process. The dynamics of M and T c over the domain are displayed in Figure 5b,e. We can see that the initial differences last for less than 200 days. At t = 200 , V = 6.25 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has higher concentration in the middle and ranges between 0.177 and 0.179 , Figure 5f. At t = 1000 , V = 34.98 V 0 (comparing with the reference case when V = 34.6 V 0 ), the distribution of C is close to uniform, Figure 5g. The initial differences in M and T c have not brought significant different behaviors in C. The steady-state value is still 1.

3.3. The Case of Using a Robin Boundary Condition for Five Types of Cells

To simulate the movements of the immune cell types: M, T c , T h , T r , and D n between the inside and outside of the tumor, we use the uniform initial value on the entire domain for all cell types and Robin boundary condition by assigning α i = 100 and α i [ X ^ i ] = 1 for these five immune cell types. Although we used the same boundary condition for all these cells, it has resulted in different distributions of T c , M, and D n over the domain at t = 200 , respectively (Figure 6a–c). The value of M over the domain varies greatly while D n becomes nearly uniform after t = 100 . At t = 200 , V = 5.09 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has lower concentration over the boundary and ranges between 0.143 and 0.147 (lower than 0.16 in the reference case), see Figure 6d. At t = 1000 , V = 33.04 V 0 (comparing with the reference case when V = 34.6 V 0 ), the dynamics of C at the center and on the boundary are displayed in Figure 6h. In this case, we observe that the cells that activate the growth of cancer cells are outperformed by cells that inhibit the growth of cancer cells. From applying the Robin boundary condition to every cell type one at a time, we observe that the changes in macrophages result in the biggest changes in the concentration of cancer cells compared to the other immune cell types.

3.4. The Case of a Positive Source of M in the Middle

We simulated the case of adding a constant source of M at all times in the middle of the tumor. The source has resulted in 200 % higher concentration of M in the middle of the tumor for most of the time interval, see Figure 7a–c. As a result, at t = 200 , V = 6.70 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has 20 % higher concentration in the middle and ranges between 0.18 and 0.21 , see Figure 7d. Moreover, at t = 1000 , V = 38.00 V 0 (comparing with the reference case when V = 34.6 V 0 ), C has a higher concentration in the middle and ranges between 1.08 and 1.15 , see Figure 7e. The dynamics of C at the center and on the boundary are displayed in Figure 7f.

3.5. The Case of a Positive Source of T c in the Middle

We also simulated the condition of a constant source of Cytotoxic and NK cells T c in the middle of the tumor. The source has resulted in a roughly 0.6 higher concentration of T c at all times, see Figure 8a,b for two time points t = 200 and t = 1000 and Figure 8c for the comparison of the maximum and minimum of T c in the whole time interval. At t = 200 , V = 4.87 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has more than 20 % higher concentration on the boundary and ranges between 0.11 and 0.14 , see Figure 8d. At t = 1000 , V = 32.30 V 0 (comparing with the reference case when V = 34.6 V 0 ), C has a lower concentration in the middle and ranges between 0.81 and 0.95 , see Figure 8e. The dynamics of C at the center and on the boundary are displayed in Figure 8f. The concentration of C is approaching its steady-state value of 1 much slower in the middle of the tumor.

3.6. The Case of a Positive Source of T c on the Boundary

We have simulated the condition of a constant source of cytotoxic and NK cells T c on the boundary of the tumor. The source has resulted in a roughly 33 % higher concentration of this cell type over the boundary at all times; see Figure 9a,b for two time points t = 200 and t = 1000 and Figure 9c for the comparison of the maximum and minimum of T c in the whole time interval. At t = 200 , V = 4.77 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has a lower concentration over the boundary and ranges between 0.127 and 0.142 (significantly lower than 0.16 in the reference case), see Figure 9d. At t = 1000 , V = 32.00 V 0 (comparing with the reference case when V = 34.6 V 0 ), C has a higher concentration in the middle and ranges between 0.89 and 0.95 , see Figure 9e. The dynamics of C at the center and on the boundary are displayed in Figure 9f.

3.7. The Case of a Positive Source of M on the Boundary

We have also simulated a constant source of macrophages M on the boundary of the tumor. The source has resulted in up to 133 % higher concentration of this cell type over the boundary at most of the time points; see Figure 10a,b for two time points t = 200 and t = 1000 and Figure 10c for the comparison of the maximum and minimum of T c in the whole time interval. At t = 200 , V = 8.88 V 0 (comparing with the reference case when V = 5.48 V 0 ), C has a lower concentration over the boundary and ranges between 0.23 and 0.27 (significantly higher than 0.16 in the reference case), see Figure 10d. At t = 1000 , V = 42.09 V 0 (comparing with the reference case when V = 34.6 V 0 ), C has a higher concentration in the middle and ranges between 1.2 and 1.25 , see Figure 10e. The dynamics of C at the center and on the boundary are displayed in Figure 10f. We can see that the cancer cell concentrations C approach different values, positively related to the concentrations of Macrophages M. Still, both M and C reach steady states.
Remark 1.
Simulations associated with all three clusters of patients’ data from the ODE paper have been done. The results are qualitatively the same. Thus, only the results obtained from using data from cluster 1 have been shown in this paper.

4. Discussion

Accumulating evidence demonstrates the critical roles of the tumor-infiltrating immune cells in tumors’ progression [28,29,30]. For example, it has been shown that cytotoxic T-cells are effector cells of adaptive immunity targeting osteosarcoma [31] and significantly affect the immune responses of osteosarcoma patients [32]. In addition, treatments with anti-tumor immunocompetence, such as NK cells and γ δ T-cells appear to be effective for some osteosarcoma tumors [33,34].
Many mathematical models have been developed to study tumors’ initiation, and progression [35,36,37,38,39,40,41,42]. Some computational models include bone modeling, osteoblast cells, or osteosarcoma treatments [43,44,45,46,47]. Recently a data-driven ODE model for the progression of osteosarcoma tumors, which considers immune cells’ interactions with tumor cells, has been developed [12]. However, to the best of our knowledge, there is no PDE model for osteosarcoma tumors considering the immune cells and tumor cell interactions. We have extended a recent ODE model for the immune and cancer cell interactions in osteosarcoma tumors [12] by developing a bio-mechanical multiphase model. Using this PDE model, which includes spatial information, namely, cellular/molecular distributions, influxes, size, and shape of the domain, we have simulated different scenarios to investigate the effect of the location of immune cells and their influx on the location and concentration of cancer cells as well as the tumor’s growth.
We have used the solution of the mechanical part, which also considers the biological variables, to trigger the change in the domain. We assess the accuracy of the domain change from the number of cells that are a part of the solid tumor. The results shown in Table 3 suggest that the domain changes fall into a reasonable range. The model also captures the killing effect of cytotoxic T-cells well. Cancer cells become more concentrated in the middle of the tumor and depleted in the boundary when there is a constant source of cytotoxic cells on the boundary of the tumor (Figure 9), and vice versa if the T c ’s are more located in the middle of the tumor, Figure 8).
The model results emphasize the importance of the influx and the location of macrophages on the cancer cells’ concentration. As we simulated the influx of immune cells, we noticed a higher level of cancer cells with the influx of macrophages than any other immune cells. Importantly, suppose there is a constant source of macrophages on the boundary of the tumors. In that case, cancer cells are collocated with them on the boundary Figure 10d, implying that the concentration of macrophages is spatially positively related to the concentration of cancer cells. This relation cannot be observed through an ODE modeling of the tumor. These results agree with the observations of a high infiltration of macrophages in localized osteosarcoma tissues using immunohistochemical staining techniques [48].
The greatest challenge in this study was finding the parameter values of the Biot equations for osteosarcoma tumors. Ideally, those parameter values should be measured exclusively for osteosarcoma. However, if we find a value for general sarcomas, we will see it as a good approximation. As such examples, we found the bulk modulus K, shear modulus G, and hydraulic conductivity κ in [19]. We found the values of the rest of the parameters from [20]; these values were calculated using an artificial neural network for general tumors. Fortunately, using those parameters, the results we obtained are biologically sound. However, the results of the model should be experimentally and bio-medically validated.

5. Conclusions

To sum up, we have proposed a model that describes osteosarcoma tumor progression and explored the effect of immune cells infiltration in the cancer cells concentration. As we are providing this model as an open-source Python code, this model can be utilized by researchers and can be improved by considering more factors, such as angiogenesis and chemotactic and haptotactic effects, or be adjusted for other cancer types. It may also be used in developing a model to elaborate the mechanical tumor growth, which could be made possible with recently published MRI image sources and image processing results [49,50].

Author Contributions

Conceptualization, L.S.; methodology, Y.H.; software, Y.H. and N.M.M.; validation, Y.H.; formal analysis, Y.H.; investigation, Y.H. and N.M.M.; resources, L.S.; data curation, Y.H.; writing—original draft preparation, Y.H., N.M.M. and L.S.; writing—review and editing, Y.H., N.M.M. and L.S.; visualization, Y.H.; supervision, L.S.; project administration, L.S.; funding acquisition, L.S. All authors have read and agreed to the published version of the manuscript.

Funding

Research reported in this publication was supported by the National Cancer Institute of the National Institutes of Health under Award Number R21CA242933. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Python scripts for computations and plotting the dynamical results are available here : https://github.com/ShahriyariLab/Bio-Mechanical-model-of-osteosarcoma-tumor-microenvironment-A-porous-media-approach, since November 2022.

Acknowledgments

The results published here are based upon the data generated by the Therapeutically Applicable Research to Generate Effective Treatments (https://ocg.cancer.gov/programs/target, accessed on 16 April 2021) initiative, phs000468.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. The ODE Model and Its Parameters

Here is the system of ODEs presented in [12] that models the biochemical interaction of key players of the osteosarcoma tumors given in Figure 1. In this paper, we have denoted the right hand-sides of this system by f i ’s. In the following equations, the proliferation and activation rates are denoted by λ ’s, and the degradation and death rates are denoted by δ ’s. In addition, the constant production rates of the cells are denoted by A’s. For the detail of the derivations of this system of ODEs, please see [12].
d M N d t = A M N λ M I γ I γ + λ M μ 1 μ 1 M N δ M N M N ,
d M d t = λ M I γ I γ + λ M μ 1 μ 1 M N δ M M ,
d T N d t = A T N λ T h M M + λ T h D D T N λ T r μ 1 μ 1 T N λ T c T h T h + λ T c M M + λ T c D D T N δ T N T N ,
d T h d t = λ T h M M + λ T h D D T N δ T h T r T r + δ T h μ 1 μ 1 + δ T h T h ,
d T r d t = λ T r μ 1 μ 1 T N δ T r T r ,
d T c d t = λ T c T h T h + λ T c M M + λ T c D D T N δ T c T r T r + δ T c μ 1 μ 1 + δ T c T c ,
d D N d t = A D N λ D C C + λ D H H D N δ D N D N ,
d D d t = λ D C C + λ D H H D N δ D C C + δ D D ,
d C d t = λ C + λ C μ 1 μ 1 + λ C μ 2 μ 2 C 1 C C 0 δ C T c T c + δ C I γ I γ + δ C C ,
d N d t = α N C δ C T c T c + δ C I γ I γ + δ C C δ N N ,
d I γ d t = λ I γ T h T h + λ I γ T c T c δ I γ I γ ,
d μ 1 d t = λ μ 1 T h T h + λ μ 1 M M + λ μ 1 C C δ μ 1 μ 1 ,
d μ 2 d t = λ μ 2 T h T h + λ μ 2 M M + λ μ 2 C C δ μ 2 μ 2 ,
d H d t = λ H M M + λ H D D + λ H N N δ H H .
Table A1. Non-dimensional initial conditions for each cluster.
Table A1. Non-dimensional initial conditions for each cluster.
Cluster M N / M N M / M T N / T N T h / T h T r / T r T c / T c D N / D N
12.3671.0050.0190.7940.7640.8281.122
20.9540.7531.2991.4512.3130.0620.071
30.8661.1040.5720.3400.48401.643
D / D C / C N / N I γ / I γ μ 1 / μ 1 μ 2 / μ 2 H / H
100.0200.1602.3941.1041.8061.059
20.6930.0050.0180.8591.3073.2590.988
300.0140.00080.2761.0301.2961.284
Table A2. Non-dimensional parameter values for each cluster, the sources can be found in the ODE paper [12].
Table A2. Non-dimensional parameter values for each cluster, the sources can be found in the ODE paper [12].
ParameterCluster 1Cluster 2Cluster 3
λ ¯ M I γ 4.3649  × 10 3 8.4234  × 10 4 2.8083  × 10 3
λ ¯ M μ 1 1.0635  × 10 2 1.4158  × 10 2 1.2192  × 10 2
λ ¯ T h M 3.3434  × 10 2 1.9270  × 10 2 2.2194  × 10 2
λ ¯ T h D 1.0963  × 10 1 7.37789.8325
λ ¯ T r μ 1 6.3  × 10 2 6.3  × 10 2 6.3  × 10 2
λ ¯ T c T h 6.11712.38462.8415
λ ¯ T c M 1.74781.22631.4683
λ ¯ T c D 1.1463  × 10 1 9.39001.3011  × 10 1
λ ¯ D C 4.0114  × 10 1 4.8942  × 10 1 5.9472  × 10 1
λ ¯ D H 4.1518  × 10 1 4.2621  × 10 1 4.1729  × 10 1
λ C 1.662  × 10 2 1.662  × 10 2 1.662  × 10 2
λ ¯ C μ 1 3.7101  × 10 4 3.5910  × 10 4 4.0692  × 10 4
λ ¯ C μ 2 7.1405  × 10 4 6.3207  × 10 4 5.7910  × 10 4
λ ¯ I γ T h 6.30951.1946  × 10 1 4.1848
λ ¯ I γ T c 2.6961  × 10 1 2.1324  × 10 1 2.9085  × 10 1
λ ¯ μ 1 T h 1.8813  × 10 2 1.1491  × 10 2 1.0315  × 10 2
λ ¯ μ 1 M 1.0751  × 10 2 1.1818  × 10 2 1.0661  × 10 2
λ ¯ μ 1 C 1.9184  × 10 2 2.5440  × 10 2 2.7772  × 10 2
λ ¯ μ 2 T h 1.23136.8806  × 10 1 6.0936  × 10 1
λ ¯ μ 2 M 1.40731.41531.2595
λ ¯ μ 2 C 2.51133.04673.2811
λ ¯ H M 4.40465.83551.8254
λ ¯ H D 3.61075.58562.0218
λ ¯ H N 5.0685  × 10 1 4.7279  × 10 1 5.4853  × 10 1
δ M N 6.93  × 10 1 6.93  × 10 1 6.93  × 10 1
δ M 1.5  × 10 2 1.5  × 10 2 1.5  × 10 2
δ T N 4.2  × 10 4 4.2  × 10 4 4.2  × 10 4
δ ¯ T h T r 6.64043.17325.0991
δ ¯ T h μ 1 4.12533.99294.5246
δ T h 2.31  × 10 1 2.31  × 10 1 2.31  × 10 1
δ T r 6.3  × 10 2 6.3  × 10 2 6.3  × 10 2
δ ¯ T c T r 1.1671  × 10 1 5.57718.9620
δ ¯ T c μ 1 7.25057.01797.9524
δ T c 4.06  × 10 1 4.06  × 10 1 4.06  × 10 1
δ D N 1.6641.6641.664
δ ¯ D C 5.3932  × 10 1 6.3864  × 10 1 7.3501  × 10 1
δ D 2.77  × 10 1 2.77  × 10 1 2.77  × 10 1
δ ¯ C T c 1.2269  × 10 2 9.6574  × 10 3 8.4017  × 10 3
δ ¯ C I γ 4.5923  × 10 3 6.4192  × 10 3 8.4660  × 10 3
δ C 3.0078  × 10 4 1.0390  × 10 3 2.4530  × 10 4
δ N 4.5935  × 10 1 4.8360  × 10 1 1.1137  × 10 1
δ I γ 3.327  × 10 1 3.327  × 10 1 3.327  × 10 1
δ μ 1 4.8748  × 10 2 4.8748  × 10 2 4.8748  × 10 2
δ μ 2 5.155.155.15
δ H 5.87  × 10 1 5.87  × 10 1 5.87  × 10 1
A ¯ M N 7.4055  × 10 1 7.0151  × 10 1 7.1382  × 10 1
A ¯ T N 1.0581  × 10 2 1.79173.1561
A ¯ D N 3.33252.39582.4867
α M N M 3.17015.6721  × 10 1 1.3878
α T N T h 1.43961.8848  × 10 1 8.8053  × 10 2
α T N T r 7.4588  × 10 1 8.2864  × 10 2 1.0267  × 10 1
α T N T c 4.65313.0144  × 10 2 1.3172  × 10 1
α D N D 2.04407.9922  × 10 1 8.1299  × 10 1

Appendix B. Details about Discretization of the PDEs

We consider a partition T h of Ω t into triangular elements and let { P j } j = 1 N h be the set of all the interior nodes of the triangulation (where N h is the number of interior nodes). Let P r denote the standard space of polynomials of total degree less than or equal to r. We introduce the space of finite elements
X h r = { g C 0 ( Ω ¯ ) : g | τ P r , τ T h } ,
that is the space of globally continuous functions that are polynomials of degree r on the single triangles of the triangulation T h . We also define
Q h = { g X h 1 : g | Ω = 0 } .
The spaces X h r and Q h are suitable for the approximation of H 1 ( Ω ) and H 0 1 ( Ω ) , respectively [51]. The space V h is selected to be ( X h 1 ) 3 . Denote the time step size by k, that is, k = T / N for some positive integer N, and t n = n k for n = 0 , 1 , , N . Let u n = u ( t n ) and p n = p ( t n ) and let U n and P n be their approximations, respectively. The fully discrete approximation for the mechanical problem is: find U n V h and P n Q h for n = 1 , 2 , , N and Lagrange multipliers λ i for i = 1 , 2 , 3 such that
a U n , v h + P n , v h + K V * n , v h i = 1 3 λ i v h , z i i = 1 3 ω i u h , z i = 0 ,
k 1 c 0 ( P n P n 1 ) , q h + k 1 α · ( U n U n 1 ) , q h + κ P n , q h = 0 ,
for all basis function v h and q h of the spaces V h and Q h . Then we solve the biological problem. Enriching the usual piece-wise linear function spaces using the bubble elements is a simple way of reducing the instability [52], for the diffusion-advection problems in the advection-dominated case [53], without increasing the size of the problem significantly [54,55].
We define
B ( l ) = span { b τ l , τ T h } ,
where b τ l are the l-th bubble functions. In this paper, we used the case r = 1 and l = 3 , thus we define the bubble enriched piece-wise linear function spaces as
B = X h 1 B ( 3 ) .
Therefore, the finite element space of the biological problem is
X h = ( B ) 14 .
The fully discrete approximation for the biological problem is: find X n ( H 1 ( Ω ) ) 14 for n = 1 , 2 , , N , such that
k 1 [ X ] n [ X ] n 1 , Ξ h + D [ X ] n , Ξ h + b [ X ] n κ Δ p , Ξ h = ( f , Ξ h ) ,
for all basis functions Ξ h of the space X h . For the specific form of cubic bubble functions or a detailed introduction of bubble elements, the readers can see Chapter 9 of [56].

Appendix C. Reports of Different Experiments

The following two types of profiles frequently show up when we plot the solution of the equations. If the variable has a higher value in the middle and a lower value over the boundary, because of different settings, as illustrated in Figure A1a, we call the profile of type I, or profile I. If, on the contrary, the variable has a higher value over the boundary and a lower value in the middle, as illustrated in Figure A1b, we call the profile of type II or profile II. Figure A1c,d show the three different locations on the domain we investigated and the dynamics of a variable at the three locations.
Figure A1. Common types of profiles in this paper. In part (a), the center has the highest density, the density reduces in the positive radial direction, and near the boundary takes its lowest value. This is because of the initial conditions of an influx through the boundary. Part (b) shows the opposite case. In part (c), three points are selected with different distance from the center. In part (d), the dynamics are plotted to show how much the density in the center is different from the boundary.
Figure A1. Common types of profiles in this paper. In part (a), the center has the highest density, the density reduces in the positive radial direction, and near the boundary takes its lowest value. This is because of the initial conditions of an influx through the boundary. Part (b) shows the opposite case. In part (c), three points are selected with different distance from the center. In part (d), the dynamics are plotted to show how much the density in the center is different from the boundary.
Cancers 14 06143 g0a1
Table A3. A report of different experiments using no-flux boundary condition.
Table A3. A report of different experiments using no-flux boundary condition.
Initial ConditionDifference Level of IC (max/min)C at t = 200 , (Relative Difference in C)V, (Relative Difference in Relevant X i ), and Others
Ref. Case, (case Section 3.1) Uniform1 C = 0.16 V = 5.48 V 0
More T h around the boundary13(< 1 % ) T h profile evened out within t = 2.5
More T c around the boundary11(< 1 % ) T c profile evened out within t = 2.5
More T c in the middle11(< 1 % ) T c profile evened out within t = 2.5
More M and T c in the middle, (case Section 3.2)11, both C [ 0.17 , 0.19 ] , profile I (Figure A1a) V = 6.25 V 0 , M maintained profile I, m a x ( T c ) / m i n ( T c ) = 2 at t = 2.5 , profile changed from I to II, then evened out
More M in the middle11Profile I V = 6.26 V 0 , M maintained profile I
More M n on the boundary30( < 1 % ) V = 5.68 V 0 , M n profile evened out at about t = 10
The relative difference is (maxmin)/min.
Table A4. A report of different experiments using uniform initial condition.
Table A4. A report of different experiments using uniform initial condition.
Boundary Condition ( α i = 1 , α i [ X ^ i ] = 1 by Default) α i , α i [ X ^ i ] C at t = 200 , (Relative Difference in C)V, (Relative Difference in Relevant X i ), and Others
Robin for all 5 cell types  --(< 0.1 % ), for all X i
Robin for all α i [ X ^ i ] = 10 maintained, ( 15 % ) V = 5.53 V 0 , M maintained profile I
Robin for T h α 4 [ X ^ 4 ] = 10 ( 0.1 % )-
Robin for T c α 6 [ X ^ 6 ] = 1000 C = 0.16 V = 5.46 V 0 , T c is 1 % more on the boundary
Robin for M--(< 1 % )
Robin for M α 2 [ X ^ 2 ] = 0.5 -(< 1 % )
Robin for M α 2 [ X ^ 2 ] = 1000 Profile II (Figure A1b), (> 10 % ) V = 7.89 V 0 , M is 200 % more on the boundary
Robin for T c --(< 0.1 % )
Robin for all α i = 10 C = 0.1543 , (< 1 % ) V = 5.45 V 0 , (≤ 0.1 % )
Robin for all,       (case Section 3.3) α i = 100 , α i [ X ^ i ] = 1 C [ 0.142 , 0.147 ] , Profile I V = 5.09 V 0 , ( [ 5 % , 15 % ] , for all i)
Robin for T c α 6 = 100 -(< 1 % )
Robin for T h α 4 = 100 -(< 1 % )
Robin for M α 2 = 100 0.151 , (< 1 % ) V = 5.09 V 0 , M [ 0.71 , 0.82 ] , profile I
The five cell types are M, Th, Tr , Tc, and Dn. Here i = 2, 4, 5, 6, 7.
Table A5. A report of different experiments using uniform initial condition and no-flux boundary condition.
Table A5. A report of different experiments using uniform initial condition and no-flux boundary condition.
SourceLocation, Max MagnitudeC at t = 200 , (Relative Difference in C)V, (Relative Difference in Relevant X i ), and Others
Gaussian for M, (case Section 3.4)In the middle, 0.1 Profile I, ( 20 % ) V = 6.70 V 0 , ( 300 % )
Gaussian for CIn the middle, 0.1 -Nonlinear solver collapsed at t = 7.5
Gaussian for CIn the middle, 0.01 -Nonlinear solver collapsed at t = 27.5
Gaussian for T c , (case Section 3.5)In the middle, 10 C [ 1.1 , 1.4 ] , profile II V = 4.87 V 0 , ( 20 % ), profile I
Gaussian for M n In the middle, 4 -Nonlinear solver collapsed at t = 70
Step function for T c On the boundary, 1Profile I, (< 5 % ) V = 5.29 V 0 , T c [ 1.3 , 1.4 ] , more on the boundary
Step function for T c , (case Section 3.6)On the boundary, 4Profile I, ( 10 % ) V = 4.75 V 0 , T c [ 1.3 , 1.7 ] , more on the boundary
Step function for T h On the boundary, 1 C = 1.573 , (< 1 % ) V = 5.57 V 0 , T h is more on the boundary, ( 8 % )
Step function for MOn the boundary, 1 C [ 0.88 , 1.1 ] , profile II V = 34.4 V 0 , M [ 5.5 , 25 ] , more on the boundary
Step function for M,       (case Section 3.7)On the boundary, 0.1 C [ 0.23 , 0.27 ] , profile II V = 8.88 V 0 , M [ 1.3 , 3.3 ] , more on the boundary

References

  1. Johns Hopkins Medicine. Osteosarcoma. Available online: https://www.hopkinsmedicine.org/health/conditions-and-diseases/sarcoma/osteosarcoma (accessed on 12 September 2020).
  2. PDQ Pediatric Treatment Editorial Board. Osteosarcoma and Undifferentiated Pleomorphic Sarcoma of Bone Treatment (PDQ®): Health Professional Version. In PDQ Cancer Information Summaries [Internet]; National Cancer Institute (US): Bethesda, MD, USA, 2022. Available online: https://www.ncbi.nlm.nih.gov/books/NBK65736/ (accessed on 8 June 2022).
  3. Hao, W.; Hauenstein, J.D.; Hu, B.; McCoy, T.; Sommese, A.J. Computing steady-state solutions for a free boundary problem modeling tumor growth by Stokes equation. J. Comput. Appl. Math. 2013, 237, 326–334. [Google Scholar] [CrossRef] [Green Version]
  4. Franks, S.; Byrne, H.; Underwood, J.; Lewis, C. Biological inferences from a mathematical model of comedo ductal carcinoma in situ of the breast. J. Theor. Biol. 2005, 232, 523–543. [Google Scholar] [CrossRef] [PubMed]
  5. Mohammad Mirzaei, N.; Tatarova, Z.; Hao, W.; Changizi, N.; Asadpoure, A.; Zervantonakis, I.K.; Hu, Y.; Chang, Y.H.; Shahriyari, L. A PDE Model of Breast Tumor Progression in MMTV-PyMT Mice. J. Pers. Med. 2022, 12, 807. [Google Scholar] [CrossRef]
  6. Shelton, S.E. Mechanistic Modeling of Cancer Tumor Growth Using a Porous Media Approach; University of North Carolina at Chapel Hill: Chapel Hill, NC, USA, 2011. [Google Scholar]
  7. Kremheller, J.; Vuong, A.T.; Yoshihara, L.; Wall, W.A.; Schrefler, B.A. A monolithic multiphase porous medium framework for (a-) vascular tumor growth. Comput. Methods Appl. Mech. Eng. 2018, 340, 657–683. [Google Scholar] [CrossRef] [PubMed]
  8. Chapman, S.J.; Shipley, R.J.; Jawad, R. Multiscale modeling of fluid transport in tumors. Bull. Math. Biol. 2008, 70, 2334–2357. [Google Scholar] [CrossRef] [PubMed]
  9. Sciume, G.; Shelton, S.; Gray, W.G.; Miller, C.T.; Hussain, F.; Ferrari, M.; Decuzzi, P.; Schrefler, B. A multiphase model for three-dimensional tumor growth. New J. Phys. 2013, 15, 015005. [Google Scholar] [CrossRef] [PubMed]
  10. Mascheroni, P.; Boso, D.; Preziosi, L.; Schrefler, B.A. Evaluating the influence of mechanical stress on anticancer treatments through a multiphase porous media model. J. Theor. Biol. 2017, 421, 179–188. [Google Scholar] [CrossRef] [Green Version]
  11. Shrestha, S.; Gokul, K.; Gurung, D.B. Temperature Variation in Breast Tissue Model with and without Tumor Based on Porous Media. J. Nepal Math. Soc. 2021, 4, 61–75. [Google Scholar] [CrossRef]
  12. Le, T.; Su, S.; Kirshtein, A.; Shahriyari, L. Data-Driven Mathematical Model of Osteosarcoma. Cancers 2021, 13, 2367. [Google Scholar] [CrossRef]
  13. Le, T.; Su, S.; Shahriyari, L. Investigating Optimal Chemotherapy Options for Osteosarcoma Patients through a Mathematical Model. Cells 2021, 10, 2009. [Google Scholar] [CrossRef]
  14. Aguda, B.; Chaplain, M.; Friedman, A.; Kimmel, M.; Levine, H.; Lolas, G.; Matzavinos, A.; Nilsen-Hamilton, M.; Swierniak, A. Tutorials in Mathematical Biosciences III: Cell Cycle, Proliferation, and Cancer; Springer: Berlin/Heidelberg, Germany, 2005. [Google Scholar]
  15. Gholami, M.D.; Saeedi, Y.; Heydari, S.; Garssen, J.; Falak, R. Exhaustion of T lymphocytes in the tumor microenvironment: Significance and effective mechanisms. Cell. Immunol. 2017, 322, 1–14. [Google Scholar] [CrossRef] [PubMed]
  16. Punt, J. Adaptive Immunity: T Cells and Cytokines. In Cancer Immunotherapy; Elsevier: Amsterdam, The Netherlands, 2013; pp. 41–53. [Google Scholar]
  17. Qiu, S.Q.; Waaijer, S.J.; Zwager, M.C.; de Vries, E.G.; van der Vegt, B.; Schröder, C.P. Tumor-associated macrophages in breast cancer: Innocent bystander or important player? Cancer Treat. Rev. 2018, 70, 178–189. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Dhodapkar, M.V.; Dhodapkar, K.M.; Palucka, A.K. Interactions of tumor cells with dendritic cells: Balancing immunity and tolerance. Cell Death Differ. 2008, 15, 39–50. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  19. Roose, T.; Netti, P.A.; Munn, L.L.; Boucher, Y.; Jain, R.K. Solid stress generated by spheroid growth estimated using a linear poroelasticity model. Microvasc. Res. 2003, 66, 204–212. [Google Scholar] [CrossRef]
  20. Dehghani, H.; Zilian, A. Poroelastic model parameter identification using artificial neural networks: On the effects of heterogeneous porosity and solid matrix Poisson ratio. Comput. Mech. 2020, 66, 625–649. [Google Scholar] [CrossRef]
  21. Hao, W.; Friedman, A. Serum upar as biomarker in breast cancer recurrence: A mathematical model. PLoS ONE 2016, 11, e0153508. [Google Scholar] [CrossRef] [Green Version]
  22. Liao, K.L.; Bai, X.F.; Friedman, A. The role of CD200–CD200R in tumor immune evasion. J. Theor. Biol. 2013, 328, 65–76. [Google Scholar] [CrossRef]
  23. Hao, W.; Friedman, A. Mathematical model on Alzheimer’s disease. BMC Syst. Biol. 2016, 10, 1–18. [Google Scholar] [CrossRef] [Green Version]
  24. Friedman, A. A free boundary problem for a coupled system of elliptic, hyperbolic, and Stokes equations modeling tumor growth. Interfaces Free Boundaries 2006, 8, 247–261. [Google Scholar] [CrossRef] [Green Version]
  25. Le, T.; Su, S.; Shahriyari, L. Immune Classification of Osteosarcoma. Math. Biosci. Eng. 2021, 18, 1879–11897. [Google Scholar] [CrossRef]
  26. Goldman, M.J.; Craft, B.; Hastie, M.; Repečka, K.; McDade, F.; Kamath, A.; Banerjee, A.; Luo, Y.; Rogers, D.; Brooks, A.N.; et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 2020, 38, 675–678. [Google Scholar] [CrossRef] [PubMed]
  27. Meyskens, F.L., Jr.; Thomson, S.P.; Moon, T.E. Quantitation of the number of cells within tumor colonies in semisolid medium and their growth as oblate spheroids. Cancer Res. 1984, 44, 271–277. [Google Scholar] [PubMed]
  28. Grivennikov, S.I.; Greten, F.R.; Karin, M. Immunity, Inflammation, and Cancer. Cell 2010, 140, 883–899. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Kitamura, T.; Qian, B.Z.; Pollard, J.W. Immune cell promotion of metastasis. Nat. Rev. Immunol. 2015, 15, 73–86. [Google Scholar] [CrossRef] [Green Version]
  30. Swann, J.B.; Smyth, M.J. Immune surveillance of tumors. J. Clin. Investig. 2007, 117, 1137–1146. [Google Scholar] [CrossRef] [Green Version]
  31. Wang, Z.; Wang, Z.; Li, B.; Wang, S.; Chen, T.; Ye, Z. Innate immune cells: A potential and promising cell population for treating osteosarcoma. Front. Immunol. 2019, 10, 1114. [Google Scholar] [CrossRef]
  32. Tsukahara, T.; Kawaguchi, S.; Torigoe, T.; Asanuma, H.; Nakazawa, E.; Shimozawa, K.; Nabeta, Y.; Kimura, S.; Kaya, M.; Nagoya, S.; et al. Prognostic significance of HLA class I expression in osteosarcoma defined by anti-pan HLA class I monoclonal antibody, EMR8-5. Cancer Sci. 2006, 97, 1374–1380. [Google Scholar] [CrossRef]
  33. Tarek, N.; Lee, D.A. Natural killer cells for osteosarcoma. Curr. Adv. Osteosarcoma 2014, 341–353. [Google Scholar] [CrossRef]
  34. Li, Z. Potential of human γδ T cells for immunotherapy of osteosarcoma. Mol. Biol. Rep. 2013, 40, 427–437. [Google Scholar] [CrossRef]
  35. Shahriyari, L.; Komarova, N.L. Symmetric vs. asymmetric stem cell divisions: An adaptation against cancer? PLoS ONE 2013, 8, e76195. [Google Scholar] [CrossRef]
  36. Shahriyari, L.; Komarova, N.L. The role of the bi-compartmental stem cell niche in delaying cancer. Phys. Biol. 2015, 12, 055001. [Google Scholar] [CrossRef] [PubMed]
  37. Shahriyari, L.; Mahdipour-Shirayeh, A. Modeling dynamics of mutants in heterogeneous stem cell niche. Phys. Biol. 2017, 14, 016004. [Google Scholar] [CrossRef] [PubMed]
  38. Bollas, A.; Shahriyari, L. The role of backward cell migration in two-hit mutants’ production in the stem cell niche. PLoS ONE 2017, 12, 2017. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Brady, R.; Enderling, H. Mathematical Models of Cancer: When to Predict Novel Therapies, and When Not to. Bull. Math. Biol. 2019, 81, 3722–3731. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  40. Chamseddine, I.M.; Rejniak, K.A. Hybrid modeling frameworks of tumor development and treatment. Wiley Interdiscip. Rev. Syst. Biol. Med. 2020, 12, e1461. [Google Scholar] [CrossRef] [Green Version]
  41. Moreia, J.; Deutsch, A. Cellular automation models of tumor development: A critical review. Adv. Complex Syst. 2002, 5, 247–267. [Google Scholar] [CrossRef]
  42. Lowengrub, J.S.; Frieboes, H.B.; Jin, F.; Chuang, Y.L.; Li, X.; Macklin, P.; Wise, S.M.; Cristini, V. Nonlinear modelling of cancer: Bridging the gap between cells and tumours. Nonlinearity 2010, 23, R1–R91. [Google Scholar] [CrossRef] [Green Version]
  43. Ji, B.; Chen, J.; Zhen, C.; Yang, Q.; Yu, N. Mathematical modelling of the role of Endo180 network in the development of metastatic bone disease in prostate cancer. Comput. Biol. Med. 2020, 117, 103619. [Google Scholar] [CrossRef]
  44. Fernández-Cervantes, I.; Morales, M.A.; Agustín-Serrano, R.; Cardenas-García, M.; Pérez-Luna, P.V.; Arroyo-Reyes, B.L.; Maldonado-García, A. Polylactic acid/sodium alginate/hydroxyapatite composite scaffolds with trabecular tissue morphology designed by a bone remodeling model using 3D printing. J. Mater. Sci. 2019, 54, 9478–9496. [Google Scholar] [CrossRef]
  45. Burova, I.; Peticone, C.; De Silva Thompson, D.; Knowles, J.C.; Wall, I.; Shipley, R.J. A parameterised mathematical model to elucidate osteoblast cell growth in a phosphate-glass microcarrier culture. J. Tissue Eng. 2019, 10, 2041731419830264. [Google Scholar] [CrossRef]
  46. Haghiralsadat, F.; Amoabediny, G.; Naderinezhad, S.; Nazmi, K.; De Boer, J.P.; Zandieh-Doulabi, B.; Forouzanfar, T.; Helder, M.N. EphA2 Targeted Doxorubicin-Nanoliposomes for Osteosarcoma Treatment. Pharm. Res. 2017, 34, 2891–2900. [Google Scholar] [CrossRef] [PubMed]
  47. Lui, G.; Treluyer, J.M.; Fresneau, B.; Piperno-Neumann, S.; Gaspar, N.; Corradini, N.; Gentet, J.C.; Marec Berard, P.; Laurence, V.; Schneider, P.; et al. A Pharmacokinetic and Pharmacogenetic Analysis of Osteosarcoma Patients Treated With High-Dose Methotrexate: Data from the OS2006/Sarcoma-09 Trial. J. Clin. Pharmacol. 2018, 58, 1541–1549. [Google Scholar] [CrossRef] [PubMed]
  48. Huang, Q.; Liang, X.; Ren, T.; Huang, Y.; Zhang, H.; Yu, Y.; Chen, C.; Wang, W.; Niu, J.; Lou, J.; et al. The role of tumor-associated macrophages in osteosarcoma progression—Therapeutic implications. Cell. Oncol. 2021, 44, 525–539. [Google Scholar] [CrossRef]
  49. Wu, J.; Yang, S.; Gou, F.; Zhou, Z.; Xie, P.; Xu, N.; Dai, Z. Intelligent segmentation medical assistance system for mri images of osteosarcoma in developing countries. Comput. Math. Methods Med. 2022, 2022, 7703583. [Google Scholar] [CrossRef] [PubMed]
  50. Lv, B.; Liu, F.; Gou, F.; Wu, J. Multi-scale tumor localization based on priori guidance-based segmentation method for osteosarcoma MRI images. Mathematics 2022, 10, 2099. [Google Scholar] [CrossRef]
  51. Quarteroni, A.; Quarteroni, S. Numerical Models for Differential Problems; Springer: Berlin/Heidelberg, Germany, 2009; Volume 2. [Google Scholar]
  52. Arnold, D.N.; Brezzi, F.; Fortin, M. A stable finite element for the Stokes equations. Calcolo 1984, 21, 337–344. [Google Scholar] [CrossRef]
  53. Brezzi, F.; Russo, A. Choosing bubbles for advection-diffusion problems. Math. Models Methods Appl. Sci. 1994, 4, 571–587. [Google Scholar] [CrossRef]
  54. Franca, L.P.; Nesliturk, A.; Stynes, M. On the stability of residual-free bubbles for convection-diffusion problems and their approximation by a two-level finite element method. Comput. Methods Appl. Mech. Eng. 1998, 166, 35–49. [Google Scholar] [CrossRef]
  55. Sendur, A. A Comparative Study on Stabilized Finite Element Methods for the Convection-Diffusion-Reaction Problems. J. Appl. Math. 2018, 2018, 4259634. [Google Scholar] [CrossRef]
  56. Gekeler, E.W. Mathematical Methods for Mechanics: A Handbook with MATLAB Experiments; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
Figure 1. Interaction network of the tumor microenvironment in osteosarcoma. Blue arrows show activation and proliferation while inhibitions are indicated by red arrows.
Figure 1. Interaction network of the tumor microenvironment in osteosarcoma. Blue arrows show activation and proliferation while inhibitions are indicated by red arrows.
Cancers 14 06143 g001
Figure 2. Dynamics of variables C and D. Part (a) shows the value of C at time t = 400 for the first group of osteosarcoma patients provided in [12]. The agreement of the curve in (b) with the green curve (results of the first group of tumors) from (c) shows that the dynamics of C at any point in the domain match the solution from the ODE paper. Sub-figure (df) show the dynamics of D.
Figure 2. Dynamics of variables C and D. Part (a) shows the value of C at time t = 400 for the first group of osteosarcoma patients provided in [12]. The agreement of the curve in (b) with the green curve (results of the first group of tumors) from (c) shows that the dynamics of C at any point in the domain match the solution from the ODE paper. Sub-figure (df) show the dynamics of D.
Cancers 14 06143 g002
Figure 3. Evolution of the domain. The original domain is plotted in (a). At t = 10 , the arrows, showing the vector field u , are pointing inward, since V is decreasing in (b) whereas in (c), they are pointing outward since V is increasing, which means that the domain is growing at t = 400 . Part (df) show the size of the domain at time t = 10 , 200, and 400, respectively.
Figure 3. Evolution of the domain. The original domain is plotted in (a). At t = 10 , the arrows, showing the vector field u , are pointing inward, since V is decreasing in (b) whereas in (c), they are pointing outward since V is increasing, which means that the domain is growing at t = 400 . Part (df) show the size of the domain at time t = 10 , 200, and 400, respectively.
Cancers 14 06143 g003
Figure 4. The diameter of the tumor is estimated from different sources. The red line, blue line, and black line represent the estimation from taking a square root, a linear model, and applying the displacement vector. It can be seen that the estimate obtained from applying the displacement vector is sound.
Figure 4. The diameter of the tumor is estimated from different sources. The red line, blue line, and black line represent the estimation from taking a square root, a linear model, and applying the displacement vector. It can be seen that the estimate obtained from applying the displacement vector is sound.
Cancers 14 06143 g004
Figure 5. Illustrations for the case of more M and T c in the middle initially. Sub-figure (a,c) shows the initial value of T c and M through the domain, respectively; (d) shows T c at t = 2.5 , which marks a drastic change in both the profile and values; the dynamics of M and T c over the whole time interval is shown in (b,e), respectively; the cancer cell concentrations at t = 100 and t = 1000 are shown in (f,g), and the dynamics of C is shown in (h).
Figure 5. Illustrations for the case of more M and T c in the middle initially. Sub-figure (a,c) shows the initial value of T c and M through the domain, respectively; (d) shows T c at t = 2.5 , which marks a drastic change in both the profile and values; the dynamics of M and T c over the whole time interval is shown in (b,e), respectively; the cancer cell concentrations at t = 100 and t = 1000 are shown in (f,g), and the dynamics of C is shown in (h).
Cancers 14 06143 g005
Figure 6. Illustrations for the case of using Robin boundary condition for 5 immune cell types. Sub-figure (ad) show the value of T c , M, D n , and C at t = 200 ; Sub-figure (eh) show the maximum and minimum in the whole interval for T c , M, D n , and C, respectively. Using the same boundary condition for 5 immune cell types results in different profiles for these cell types.
Figure 6. Illustrations for the case of using Robin boundary condition for 5 immune cell types. Sub-figure (ad) show the value of T c , M, D n , and C at t = 200 ; Sub-figure (eh) show the maximum and minimum in the whole interval for T c , M, D n , and C, respectively. Using the same boundary condition for 5 immune cell types results in different profiles for these cell types.
Cancers 14 06143 g006
Figure 7. Illustrations for the case of a source of M in the middle. Sub-figure (a,b) show the value of M at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of M, from where we can see that the source has introduced a significant amount of increase in the middle; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f). These figures suggest that the cancer cells will be at the place where there are more macrophages, yet the concentration still reaches a steady state.
Figure 7. Illustrations for the case of a source of M in the middle. Sub-figure (a,b) show the value of M at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of M, from where we can see that the source has introduced a significant amount of increase in the middle; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f). These figures suggest that the cancer cells will be at the place where there are more macrophages, yet the concentration still reaches a steady state.
Cancers 14 06143 g007
Figure 8. Illustrations for the case of a source of T c in the middle. Part (a,b) show the value of T c at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of T c , depicting that the value of T c in the middle is constantly higher than it on the boundary by a considerable amount; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f), where the growth of cancer cells in the middle are inhibited strongly by the higher concentration of T c .
Figure 8. Illustrations for the case of a source of T c in the middle. Part (a,b) show the value of T c at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of T c , depicting that the value of T c in the middle is constantly higher than it on the boundary by a considerable amount; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f), where the growth of cancer cells in the middle are inhibited strongly by the higher concentration of T c .
Cancers 14 06143 g008
Figure 9. Illustrations for the case of a constant source of T c on the boundary. Su-figures (a,b) shows the value of T c at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of T c , depicting that the value of T c on the boundary is constantly higher than it is elsewhere; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f), where the growth of cancer cells, especially near the boundary, are inhibited strongly by the higher concentration of T c .
Figure 9. Illustrations for the case of a constant source of T c on the boundary. Su-figures (a,b) shows the value of T c at t = 200 and t = 1000 , respectively; (c) shows the dynamics of the maximum and minimum of T c , depicting that the value of T c on the boundary is constantly higher than it is elsewhere; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f), where the growth of cancer cells, especially near the boundary, are inhibited strongly by the higher concentration of T c .
Cancers 14 06143 g009
Figure 10. Illustrations for the case of a source of M on the boundary. Sub-figure (a,b) show the value of T c at t = 200 and t = 1000 , respectively; (c) show the dynamics of the maximum and minimum of M, depicting that the value of M on the boundary is substantially higher than it is elsewhere; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f). We can see that the cancer cell concentrations approach different steady-state values, positively related to the concentrations of Macrophages.
Figure 10. Illustrations for the case of a source of M on the boundary. Sub-figure (a,b) show the value of T c at t = 200 and t = 1000 , respectively; (c) show the dynamics of the maximum and minimum of M, depicting that the value of M on the boundary is substantially higher than it is elsewhere; the cancer cell concentration at t = 200 and t = 1000 are shown in (d,e); the value of the maximum and minimum of C over the whole time interval is shown in (f). We can see that the cancer cell concentrations approach different steady-state values, positively related to the concentrations of Macrophages.
Cancers 14 06143 g010
Table 1. Variable names corresponding to [ X i ] .
Table 1. Variable names corresponding to [ X i ] .
[ X i ] Variable NameBiological Meaning (Concentration of)Scaling Factor
[ X 1 ] M n Naive macrophages 6.236 × 10 6
[ X 2 ] MMacrophages 1.977 × 10 7
[ X 3 ] T n Naive T cells 4.926 × 10 6
[ X 4 ] T h Helper T cells 7.092 × 10 6
[ X 5 ] T r Regulatory T cells 3.675 × 10 6
[ X 6 ] T c Cytotoxic T cells and NK cells 2.292 × 10 7
[ X 7 ] D n Naive dendritic cells 4.826 × 10 5
[ X 8 ] DDendritic cells 9.865 × 10 5
[ X 9 ] CCancer cells 1.343 × 10 10
[ X 10 ] NNecrotic cells 3.764 × 10 8
[ X 11 ] I γ IFN- γ 0.868
[ X 12 ] μ 1 TGF- β , IL-4, IL-10, and IL-13 21.510
[ X 13 ] μ 2 IL-6 and IL-17 2.067
[ X 14 ] H n HMGB1 5.076
These are values associated with cluster 1 in the ODE paper [12].
Table 2. Mechanical parameter values.
Table 2. Mechanical parameter values.
ParameterNameValue
ϕ Porosity0.2
KBulk modulus40,000 (Pa) [19]
GShear modulus30,000 (Pa) [19]
α Biot effective stress coefficient0.7 [20]
κ Hydraulic conductivity 6.9 × 10 14 ( m 2 · Pa 1 · s 1 ) [19]
MBiot modulus 2 × 10 5 (Pa) [20]
D c e l l Diffusion coefficient for cells 3.6 × 10 8 ( cm 2 · h 1 ) [21]
D c y t o Diffusion coefficient for cytokines 5.2 × 10 5 ( cm 2 · h 1 ) [22]
D H Diffusion coefficient for HMGB1 3.3 × 10 3 ( cm 2 · h 1 ) [23]
These values are for sarcoma, not necessarily osteosarcoma; These values are for generic tumors, evaluated through Artificial Neural Network (ANN) and then averaged among all reasonable values.
Table 3. The total number of cells V, the ratios n = V / V 0 , the crude estimate of diameters 0.01 n , the linear model’s prediction of diameters, and the diameter of the domain after we apply the solid displacement vector u , at different times.
Table 3. The total number of cells V, the ratios n = V / V 0 , the crude estimate of diameters 0.01 n , the linear model’s prediction of diameters, and the diameter of the domain after we apply the solid displacement vector u , at different times.
t Cell Number V  Ratio n = V / V 0 Crude Estimate 0.01 n Linear Model EstimateSimulated Domain Diameter
033,72610.010.010.01
1024,2640.7850.00890.00870.0096
5041,4471.2290.01110.01120.0102
10071,7782.1280.01460.01530.0116
150118,1723.5040.01870.02020.014
200184,5785.4730.02340.02590.0176
250273,5108.1100.02850.03240.0227
300384,15611.390.03370.03920.0291
350510,78115.140.03890.04590.0367
400642,92519.060.04370.05230.0451
450768,23322.780.04770.05780.0535
500876,63025.990.05100.06220.0614
550963,10928.560.05340.06560.0684
6001,027,71130.470.05520.06800.0743
6501,073,62231.830.05640.06970.0793
7001,105,10432.770.05720.07090.0834
7501,126,16533.390.05780.07160.0868
8001,140,02233.800.05810.07210.0897
8501,149,04234.070.05840.07240.0922
9001,154,87234.240.05850.07260.0944
9501,158,62234.350.05860.07280.0963
10001,160,21534.400.05870.07280.098
These values are cell numbers in a 2D area. A scale of 106 in a 2D area amounts to a scale of 109 in a 3D volume.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hu, Y.; Mohammad Mirzaei, N.; Shahriyari, L. Bio-Mechanical Model of Osteosarcoma Tumor Microenvironment: A Porous Media Approach. Cancers 2022, 14, 6143. https://doi.org/10.3390/cancers14246143

AMA Style

Hu Y, Mohammad Mirzaei N, Shahriyari L. Bio-Mechanical Model of Osteosarcoma Tumor Microenvironment: A Porous Media Approach. Cancers. 2022; 14(24):6143. https://doi.org/10.3390/cancers14246143

Chicago/Turabian Style

Hu, Yu, Navid Mohammad Mirzaei, and Leili Shahriyari. 2022. "Bio-Mechanical Model of Osteosarcoma Tumor Microenvironment: A Porous Media Approach" Cancers 14, no. 24: 6143. https://doi.org/10.3390/cancers14246143

APA Style

Hu, Y., Mohammad Mirzaei, N., & Shahriyari, L. (2022). Bio-Mechanical Model of Osteosarcoma Tumor Microenvironment: A Porous Media Approach. Cancers, 14(24), 6143. https://doi.org/10.3390/cancers14246143

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop