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

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.


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.

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.

Biological Part of the Model
We adopt the Reaction-Diffusion-Advection equations by adding the diffusion term D cell ∆[X i ] and the advection term b∇ · ([X i ]κ∇p) to the ODEs of each cell type. Here, D cell 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 cyto .
In summary, the biological part of our model can be written as where D i is 0, D cell , D cyto , 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 where D = diag{D i } and B = diag{b i } are diagonal coefficient matrices.

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 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.

Governing Equations
We start with presenting the stress(σ ij )-strain(e ij ) relation for linear, isotropic poroelastic material. In the spirit of Hooke's law, the linear relation is Here, e ij = 1 2 volumetric strain, where u is the solid displacement of the material, δ ij 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: The volume of cells that are a part of the solid tumor V is calculated through 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, 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: where we define u as the solid displacement of the solid tumor, a macroscopic measure.
Constitutive equation of the fluid pressure p: where M is the Biot modulus. Equilibrium equation (Newton's law of motion): Darcy's law (flow of a fluid with constant viscosity across a solid is a function of its pressure difference and the solid properties): 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 Solving the governing Equations (7)-(12) simultaneously, we obtain the mechanical part of our model: In Ω(t) which is to be aligned with the biological part (1) 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.

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 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 ) 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.
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] ∂Ω udx = 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 or, 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.

Weak Formulation
We first rewrite the biological equations in the Eulerian system. Noticing that the definition of the material derivative is given by we have where v f is the fluid velocity. Let I = (0, T], V = (L 2 (I; H 1 (Ω))) 3 , Q = L 2 (I; H 1 (Ω)), (16) and Lagrange multipliers λ i such that for a.e. t ∈ (0, T] and 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: where z i ∈ {(1, 0), (0, 1), (−x 2 , x 1 )}.

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 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 for all basis functions Ξ h of the space X h .

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 Sections 3.4-3.7, as shown below.

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: p(x, t) = p 0 , on ∂Ω, and initial conditions: 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). 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 nV 0 , we have an intuitive but crude estimate of the tumor size at that time, A = nA 0 = π 4 × 10 −4 n (m 2 ), and an estimate of the diameter of the tumor d = √ nd 0 (since we are working on 2D cases). For example, by the time the total number of cells becomes 4V 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 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.

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.25V 0 (comparing with the reference case when V = 5.48V 0 ), C has higher concentration in the middle and ranges between 0.177 and 0.179, Figure 5f. At t = 1000, V = 34.98V 0 (comparing with the reference case when V = 34.6V 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.

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.09V 0 (comparing with the reference case when V = 5.48V 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.04V 0 (comparing with the reference case when V = 34.6V 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.

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.87V 0 (comparing with the reference case when V = 5.48V 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.30V 0 (comparing with the reference case when V = 34.6V 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.

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

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.88V 0 (comparing with the reference case when V = 5.48V 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.09V 0 (comparing with the reference case when V = 34.6V 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.

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.

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].

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].

Appendix B. Details about Discretization of the PDEs
We consider a partition T h of Ω t into triangular elements and let {P j } N h j=1 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 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 The spaces X r h and Q h are suitable for the approximation of H 1 (Ω) and H 1 0 (Ω), respectively [51]. The space V h is selected to be (X 1 h ) 3 . Denote the time step size by k, that is, k = T/N for some positive integer N, and t n = nk 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 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 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 Therefore, the finite element space of the biological problem is The fully discrete approximation for the biological problem is: find X n ∈ (H 1 (Ω)) 14 for n = 1, 2, . . . , N, such that 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. Table A3. A report of different experiments using no-flux boundary condition.  Step function for T c On the boundary, 1 Profile I, (< 5%) V = 5.29V 0 , T c ∈ [1.3, 1.4], more on the boundary