A mechanobiological model for tumor spheroid evolution with application to glioblastoma: a continuum multiphysics approach

Computers in Biology and Medicine 159 (2023) 106897 Available online 18 April 2023 0010-4825/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/compbiomed A mechanobiological model for tumor spheroid evolution with application to glioblastoma: A continuum multiphysics approach Ana Carrasco-Mantisa, Teodora Randelovicb,c, Héctor Castro-Abrilb,c,d, Ignacio Ochoab,c,d, Manuel Doblaréb,c,d, José A. Sanz-Herreraa,∗ aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Spain bAragón Institute of Engineering Research (I3A), University of Zaragoza, Spain cAragón Institute of Health Research (IIS), Spain dCentro de Investigación Biomédica en Red en Bioingeniería, Biomateriales y Nanomedicina (CIBER-BBN), Spain ARTICLE INFO Keywords: Cellular spheroid Glioblastoma Mechanobiology Numerical simulation Finite element method ABSTRACT Background: Spheroids are in vitro quasi-spherical structures of cell aggregates, eventually cultured within a hydrogel matrix, that are used, among other applications, as a technological platform to investigate tumor formation and evolution. Several interesting features can be replicated using this methodology, such as cell communication mechanisms, the effect of gradients of nutrients, or the creation of realistic 3D biological structures. The main objective of this work is to link the spheroid evolution with the mechanical activity of cells, coupled with nutrient consumption and the subsequent cell dynamics. Method: We propose a continuum mechanobiological model which accounts for the most relevant phenomena that take place in tumor spheroid evolution under in vitro suspension, namely, nutrient diffusion in the spheroid, kinetics of cellular growth and death, and mechanical interactions among the cells. The model is qualitatively validated, after calibration of the model parameters, versus in vitro experiments of spheroids of different glioblastoma cell lines. Results: Our model is able to explain in a novel way quite different setups, such as spheroid growth (up to six times the initial configuration for U-87 MG cell line) or shrinking (almost half of the initial configuration for U-251 MG cell line); as the result of the mechanical interplay of cells driven by cellular evolution. Conclusions: Glioblastoma tumor spheroid evolution is driven by mechanical interactions of the cell aggregate and the dynamical evolution of the cell population. All this information can be used to further investigate mechanistic effects in the evolution of tumors and their role in cancer disease. 1. Introduction Spheroids are increasingly used in the recent years as an enabling biological platform to study tumor formation, evolution, as well as to study cancer treatments [1]. Spheroids are clusters of cells, eventually cultured within a hydrogel matrix that mimics the extracellular matrix, which show a quasi-spherical morphology in vitro [2]. Even though spheroids lack the enormous complexity of tumor tissues, they are useful in vitro models as they are able to replicate several important features of solid tumors. In particular, communication mechanisms and cell differentiation have been observed in spheroids [3]. Also, the high cell density and compaction in larger spheroids induce areas that suffer hypoxia and gradients of nutrients, as the diffusion within the structure is limited, becoming a good surrogate for these effects that happen in reality [4]. Moreover, the 3D environment of the spheroid, versus 2D, allows a better understanding of the mechanisms of cellular ∗Correspondence to: Camino de los descubrimientos s/n, 41092 Seville, Spain. E-mail address: [email protected] (J.A. Sanz-Herrera). communication and signaling, both autocrine and paracrine, and the different types of cell migration stimuli, such as chemotaxis, durotaxis, electrotaxis and haptotaxis [5]. Tumor heterogeneity (intertumoral or intratumoral) is another important characteristic, related to resistance, which can be also reproduced by spheroids [6]. From an in vitro point of view, spheroids have been used to study the mechanisms of cell differentiation, cell contact and the different cell phenotypes that give rise to different regions within the spheroid [7]. Furthermore, Friedrich et al. [8] developed platforms to study anticancer treatments in spheroids. Nath and Devi [9] presented several spheroid generation techniques and assays for their characterization. On the other hand, Mark et al. [10] studied the relationship between the cellular mechanics of the spheroid and the deformation of the extracellular matrix; and Ayuso et al. [11] established that spheroids https://doi.org/10.1016/j.compbiomed.2023.106897 Received 11 November 2022; Received in revised form 9 February 2023; Accepted 9 April 2023 Computers in Biology and Medicine 159 (2023) 106897 2 A. Carrasco-Mantis et al. represent a good tool in the study of migration by chemotaxis, the main mechanism of cell migration. In another study, Guillaume et al. [12] suggested that the stresses that originate in a spheroid during its growth affect its physical properties and therefore the response to treatment. Wu et al. [13] developed microfluidic devices for the study of chemotaxis, through the generation of chemical gradients in both simple and complex environments, and studied chemotaxis as a function of the cell type. The devices designed by Cheng et al. [14] allow a quick reconfiguration for use in different applications, including the analysis of tumor first development and evolution. Ayuso et al. [15] performed a protocol for the confinement of the gel in microfluidic chips and, subsequently [16], they were able to reproduce the complex tumor environment in these platforms. Fatehullah et al. [17] presented a generation of organoids from biopsies performed on tumor tissues and their subsequent cultivation in hydrogels, which constitute a novel test in the study of pathologies and treatments. From an in silico perspective, Xue et al. [18] developed a continuum model to investigate the relationship between biology, chemistry and mechanics in the evolution of tumor spheroids by treating the tumor as a porous system. They highlighted the importance of nutrient transport and mechanical stresses in growth. The model was validated with data taken from the literature. Also, Mascheroni et al. [19] presented a continuum tumor spheroid growth model whose equations were developed according to the porous media theory. This study also performed experiments with the U-87 MG tumor line, including both free and pressurized growth. However, the diffusion of nutrients into the different cell types of the spheroid was not considered. Greenspan [20] developed a continuum model to analyze growth of a colony of cells, which compares to the movement of an incompressible fluid in a membrane of variable size. Growth was considered to be influenced by a surface tension force and by the diffusion of nutrients in this work. Nonetheless, the analysis was not conducted for spheroids, and the results were not compared with experiments. Byrne and Chaplain [21] showed that nutrient concentration at the spheroid boundary is not constant, but follows the Gibbs–Thomson relationship, which states that the local curvature at the edge of the spheroid influences this concentration, carrying out a mathematical description of this phenomenon. The results were not validated in this study. Furthermore, Chen et al. [22] proposed a model that analyzes how the mechanics of the elastic and isotropic extracellular medium affects tumor evolution, but it does not model the mechanical interactions of the cells in the spheroid. The authors stated that when a spheroid grows, it exerts a pressure on the medium, which in turn exerts a stress on the spheroid contour. In addition, they performed an experimental validation of their model, based on results found in the literature. Ambrosi and Mollica [23] described the growth of a tumor spheroid, influenced by nutrient diffusion and biochemical factors, considering the same number of particles at each instant of time and justifying growth as an increase of mass of the particles but not on their number. Hence, the dynamics of cellular proliferation and death was not considered, and an experimental validation of the model was not conducted. In addition, Ambrosi and Mollica [24] studied the growth of a spheroidal tumor as a soft tissue and highlighted the need to know the history of inhomogeneous mass increase. The authors developed a mathematical model for the experiments performed by other authors. In a different study, Ambrosi et al. [25] assumed that the spheroid behaves as a poroelastic material whose growth is influenced by the stress applied on its boundary. The authors also studied the presence of residual stresses in the spheroid. In order to experimentally validate their model, they based their work on previous existing studies in the literature. Fraldi and Carotenuto [26] carried out a poroelastic model in which an inelastic growth derived from a Volterra/Lotka type interaction between different cell types occurs. However, nutrient diffusion was not considered in this work. Additionally, Bull et al. [27] recently proposed a discrete (agent-based) approach to model the dynamics of cells in spheroids, accounting for the diffusion of nutrients and oxygen from the medium to the interior of the spheroid, which determines the cellular actions of proliferation, differentiation and cell death. In the same context, Amereh et al. [28] presented a mathematical model of tumor formation in spheroids based on reaction–diffusion equations. With regard to mathematical mechanobiological models, there is a vast literature on different methodologies and mathematical approaches to study in silico the effect of mechanics in biological processes. To cite a few, Perfahl et al. [29] developed a vasculogenesis model which considers the application of forces through the interaction of cells with their neighbors and with the extracellular environment, obtaining a vascular network as results of the simulations. Ayensa- Jiménez et al. [30], described a mathematical model of cell dynamics, growth and death, in a microfluidic device with application to glioblastoma evolution. Laird [31] performed a mathematical model of tumor proliferation with different growth stages. Furthermore, Pettet et al. [32] implemented a mathematical model of cell migration in spheroids, where chemotaxis plays an important role. Lu and Kang [33] analyzed how hypoxia-inducing factors are related to increased treatment failure and mortality, while Curtis and Seehar [34] studied the relationship between cell mechanics and proliferation, concluding that compression increases proliferation. Moreo et al. [35] carried out a predictive mechanosensing model that explains the interaction between the cell cytoskeleton and the extracellular matrix through focal adhesions. In addition, they also studied cell migration by modeling the evolution of cell concentration from a continuum point of view. Furthermore, Sanz-Herrera et al. [36] studied how, through contractile forces, the contact between a cell and the substrate occurs when it shows a curvature, analyzing the influence of microfilament bending and the initial shape of the substrate. Huang and Ingber [37] described how through integrins, focal adhesions and Rho/ROCK signaling that enables cellular contractility, a mechanistic process of autocrine cell communication important in cancer development, occurs. Additionally, Nagelkerke et al. [38] highlighted the role of mechanics in cancer development, describing several processes, including activation of the Rho family of proteins, involved in intracellular contractility of the cytoskeleton. Our application is devoted to analyze glioblastoma spheroids. Glioblastoma is the most malignant brain tumor, with a life expectancy after the standard treatment (surgery plus radio and/or chemotherapy) of about 14 months from diagnosis [39]. It originates in the glial cells and its initial evolution is controlled by thrombosis, hypoxia, cell migration events that occur successively [40] and is characterized by microvascular hyperplasia and necrosis. The area where the latter takes place (necrotic core), is usually surrounded by a high cell density, known as ‘pseudopalisade’. The hypoxic tissue (more resistant than a healthy one [41]) increases tumor glycolysis to provide energy to the cell, promotes angiogenesis, that is the creation of blood vessels and nourishes the tumor invasion and metastasis by the action of hypoxia-inducible factors [33]. In this article, we develop a phenomenological mechanobiological model to investigate the in vitro evolution of suspended glioblastoma spheroids. Compared to previous existing models in the literature, our multiphysics approach describes the evolution of spheroids from the internal mechanical activity of cells in response to nutrient stimuli in the microenvironment in a continuous way. Therefore, the model takes into consideration the main physics observed during the process, such as nutrient diffusion within the spheroid, dynamics of cellular growth and death, and mechanical interactions of cells including computation of biomechanical forces. These latter forces present at the cellular level may have an external origin or may be generated by the cells themselves. They play a fundamental role in maintaining homeostasis in the organism. They can also contribute to the development of pathologies, such as tumor progression, by increasing the rigidity of the extracellular matrix or modifying the flow of interstitial fluid in the matrix [42]. In this work, we test how biomechanical stresses and deformations, generated at the cellular level, influence the organization and morphology Computers in Biology and Medicine 159 (2023) 106897 3 A. Carrasco-Mantis et al. Fig. 1. Multiphysics model: Schematics of the different physics considered in the modeling of suspended spheroid evolution. Source: Figure created with BioRender.com. Table 1 Multiphysics mechanobiological model quantities. Symbol Description Units 𝑐(x, 𝑡)Nutrient concentration in the spheroid domain 𝑚𝑜𝑙 𝑚𝑚3 𝑛(x, 𝑡)Cell concentration in the spheroid domain (superscript g-growth, d-death) 𝑐𝑒𝑙𝑙 𝑚𝑚3 𝜀𝑇(x, 𝑡)Total strain at a material point of the spheroid – 𝜀𝑒𝑙 (x, 𝑡)Elastic strain at a material point of the spheroid (passive component) – 𝜀𝑐𝑒𝑙𝑙 (x, 𝑡)Cell strain at a material point of the spheroid as a consequence of internal cell contractility (active component) – 𝜎(x, 𝑡)Stress at a material point of the spheroid as a consequence of cell contractility 𝑃 𝑎 Table 2 Model parameters. Symbol Description Units 𝑐0Nutrient concentration in the spheroid domain at time 0 mol mm3 𝑐𝛤Nutrient concentration on the spheroid boundary mol mm3 𝑘𝑐Coefficient of nutrient diffusion from the medium to the spheroid mm2 min 𝑟Nutrient consumption rate mol cell⋅min 𝑛0Cell concentration in the spheroid domain at time 0 cell mm3 𝑐𝑛Nutrient concentration threshold for cells to survive and grow mol mm3 𝑐𝑠𝑎𝑡 Nutrient concentration saturation value which defines the maximum rate level achieves by cells to divide and grow mol mm3 𝑘𝑑Cell death rate mm3 mol⋅min 𝑘𝑔Cell growth rate mm3 mol⋅min 𝑛∗Cell concentration from which cell contractility changes from compression to expansion cell mm3 𝑘𝜖Cell expansion coefficient as consequence of cell forces exerted by the cytoskeleton machinery mm3 cell 𝛼Shape factor coefficient of nutrients in the suspension experiment mol cell⋅min of the tumor spheroids. The overall continuum model was discretized in time and is formulated in the updated Lagrangian configuration for each time step, being spatially discretized by means of the finite element method. The model permits to predict the nonlinear evolution of the spheroid along time, in terms of a number of model parameters. This dependence is investigated in a parametric analysis of the model. Computers in Biology and Medicine 159 (2023) 106897 4 A. Carrasco-Mantis et al. Fig. 2. Cellular growth/death kinetics represented as a function of the nutrient concentration. Fig. 3. Active cell’s contractility/expansion deformation model as a function of the cell concentration, for death (D) and growth (G) fates. 𝑛∗is a model parameter that represents the concentration level for which a new daughter cell cannot be accommodated within the reference volume. 𝑛0represents the initial concentration of cells. Finally, two experimental tests were designed ad-hoc in this study to qualitatively validate the proposed model. Specifically, the evolution of 1000-cell spheroids of glioblastoma cell lines U-251 MG and U-87 MG were performed and recorded. The model fits the observed trends in the experiments for a set of calibrated parameters, suggesting that the implemented physics can explain the observed behavior. In particular, it is concluded that spheroid growth or shrinking is the result of the mechanical interplay of cells driven by cellular evolution. The paper is organized as follows: first the experimental setup and the mathematical model are introduced in Sections 2–4. The numerical implementation of the model is elaborated in Section 5, and the parametric analysis is presented in Section 6.1. The validation of the proposed model versus the obtained experimental results is described in Section 6.2. Finally, the discussion of the results and main derived conclusions are drawn at the end of the paper. 2. Experimental setup Glioblastoma cell lines, U-251 MG and U-87 MG, were purchased from Sigma Aldrich and cultured in high glucose Dulbecco’s modified Eagle’s medium (DMEM; Lonza BE12-614F), supplemented with 10% fetal bovine serum (FBS; Sigma Aldrich F7524), 2 mM L-glutamine (Lonza 17-605C) and penicillin/streptomycin (Lonza 17-602E). Cell cultures were maintained at 37 ◦C within a humidified TEB-1000 incubator (EBERS Medical Technology) with 5% CO2. Spheroids were generated by hanging drop method. Cells were trypsinized and resuspended in growth medium supplemented with 20% methylcellulose solution to reach the concentration of 40 000 cells/mL. In order to form 1000-cell spheroids, 25 uL drops were placed on the lid of a Petri dish. The bottom part of the dish was filled with distilled water to prevent evaporation. Plates were placed within the incubator and left for 48 h to assure spheroid formation. Afterwards, spheroids were transferred to growth medium in suspension 96 well plates (Sarstedt 83.3925.500) treated with anti-adherence solution (Stemcell 07010). Their behavior was followed by phase contrast microscopy (Nikon TiE) during two weeks. Images were analyzed with a Fiji plugin SpheroidJ, which allows automated quantification of the spheroid area [43]. Growth curves were obtained by normalizing those values to the area that spheroids occupied at the beginning of the experiment (day 0). 3. Mathematical model The proposed mathematical model considers the main phenomena that take place in suspended spheroid evolution, according to experimental observations. In particular, spheroid size evolution is considered to be mainly determined by the proliferation of cells as well as the contractile activity exerted by cells [11,44]. Moreover, nutrient diffusion within the spheroid enhances the dynamics of cells [44], whereas lack of nutrients promotes the formation of a necrotic core in the spheroid [11]. According to the referred evidences the physics assumed by the model are fundamentally three. First, diffusion of soluble substances in the medium, such as nutrients and oxygen, to the spheroid interior. Second, cellular activity and dynamics in terms of proliferation and death, which in turn determine spheroid growth or shrinking evolution from the initial state. Finally, we model cellular contractility and mechanical interactions of cell organization in the spheroid. A spherical geometry of the spheroid is assumed, similar to other works [18,25,26]. Furthermore, as main assumptions, we consider that (i) cell dynamics depends on nutrient concentration and (ii) volume variation of the spheroid is related to cell contractility and concentration. The overall model turns into a mechanobiological and multiphysics coupled approach which is presented in a dimensionless fashion in Section 4. The schematics of the different physics of the model can be seen in Fig. 1. Model quantities and parameters are shown in Tables 1 and 2, respectively. The model is also elaborated in a numerical framework in Section 5. 3.1. Nutrient diffusion We consider that soluble substances in the medium diffuse towards the spheroid core, leading to nutrient supply to interior cells (see Fig. 1A). Therefore, nutrients (considering nutrients as oxygen and glucose concentration) are assumed, as a first approach, to follow aFickean diffusion model, analogously to other works in a similar context [34,45–47]: 𝜕𝑐 𝜕𝑡 =𝑘𝑐▽2𝑐−𝑟⋅𝑛in 𝛺(𝐱, 𝑡) 𝑐(𝐱,0) = 𝑐0in 𝛺(𝐱,0) 𝑐(𝐱, 𝑡) = 𝑓(𝐱, 𝑡)in 𝛤𝑏(𝐱, 𝑡) (1) The right-hand side of Eq. (1)a also includes a second (source) term which accounts for the consumption rate of nutrients by cells. 𝑘𝑐 is the diffusion coefficient of the spheroid aggregation, 𝑟is the rate of nutrient consumption by cells, and 𝑐0the initial concentration of nutrients within the spheroid. 𝛺represents the domain and 𝛤𝑏the contour of the spheroid. On the other hand, the function 𝑓represents Computers in Biology and Medicine 159 (2023) 106897 5 A. Carrasco-Mantis et al. Table 3 Dimensionless values of the model parameters for the parametric analysis. Parameter Reference value Lowest value Highest value  𝑘𝑐0.01 0.01 0.01 𝑟 1⋅10−2 5⋅10−3 2⋅10−2  𝑘𝑑1.5⋅10−4 0.75 ⋅10−4 3⋅10−4  𝑘𝑔1.5⋅10−4 0.75 ⋅10−4 3⋅10−4  𝑘𝑒𝑑 4⋅10−2 2⋅10−2 8⋅10−2  𝑘𝑒𝑔 4⋅10−2 2⋅10−2 8⋅10−2 𝑛∗2 1 4 𝑐𝑛0.5 0.5 0.5 𝑐𝑠𝑎𝑡 0.9 0.9 0.9 the consumption of nutrients available in the suspension experiment, defined as: 𝜕𝑓 𝜕𝑡 = −𝛼⋅𝑛in 𝛤𝑏(𝐱, 𝑡) 𝑓(𝐱,0) = 𝑐𝛤in 𝛤𝑏(𝐱,0) (2) with 𝛼a factor (shape) coefficient which accounts for the amount and availability of nutrients and configuration of the suspension experiment for diffusion of nutrients. This coefficient is introduced as a correction factor to the full and plenty availability of nutrients hypothesis in the experiment. 𝑐𝛤represents the nutrient concentration in the suspension at the initial time of the analysis. Note that, according to Eq. (2), in our model, nutrient concentration may become very low (even negative) in the medium for large time 𝑡. Indeed, this behavior is aligned with experimental observations where a recharging/refreshing protocol is needed before complete depletion of nutrients within the medium. This experimental issue can be properly introduced in our model by assuming an update of the function 𝑓at certain time points. Of course, negative concentrations do not have any biophysical meaning so our solution would not be valid under this situation, as the complete depletion of nutrients would have been achieved before. In a different way, Mascheroni et al. [19], alternatively to a Fickean diffusion as considered in Eq. (1), established a mass exchange for nutrient consumption by tumor cells; and Ambrosi et al. [25] compared spheroid diameter with nutrient penetration length to differentiate between volumetric or surface growth. 3.2. Cell dynamics Spheroid growth or decrease is ultimately depends on the dynamics of the cellular activity (see Fig. 1B). Therefore, the dynamics of cellular growth/death is represented, similar to other works [31,54,55], by means of the following first order differential equation: 𝜕𝑛 𝜕𝑡 =𝑘𝑛(𝑐)⋅𝑛in 𝛺(𝐱, 𝑡) 𝑛(𝐱,0) = 𝑛0in 𝛺(𝐱,0) (3) In Eq. (3), and as a simplification of other models [30,33,56], the kinetics of cellular growth/death is characterized by the function 𝑘𝑛(𝑐), which is assumed to be dependent of the nutrient concentration [30]. 𝑘𝑛(𝑐)is the function shown in Fig. 2, where two regions can be identified: (i) a cell death (or dormancy) region which appears for low nutrient concentration, and (ii) a cell growth zone for high levels of nutrients. Both regions are regulated by a threshold parameter 𝑐𝑛, which represents the homeostatic level of nutrients in cell’s activity. We consider a linear model for 𝑘𝑛for both death and growth regions, finishing this latter with a constant region after a certain saturation level of nutrients 𝑐𝑠𝑎𝑡. Therefore, Eq. (3) can be split into two different conditions for death and growth dynamics as follows: 𝜕𝑛𝑑 𝜕𝑡 = −𝑘𝑑⋅(𝑐𝑛−𝑐)⋅𝑛if 𝑐≤𝑐𝑛in 𝛺(𝐱, 𝑡) 𝑛𝑑(𝐱,0) = 0 in 𝛺(𝐱,0) (4) 𝜕𝑛𝑔 𝜕𝑡 =𝑘𝑔⋅(𝑐−𝑐𝑛)⋅𝑛if 𝑐 > 𝑐𝑛in 𝛺(𝐱, 𝑡) 𝜕𝑛𝑔 𝜕𝑡 =𝑘𝑔⋅(𝑐𝑠𝑎𝑡 −𝑐𝑛)⋅𝑛if 𝑐 > 𝑐𝑠𝑎𝑡 in 𝛺(𝐱, 𝑡) 𝑛𝑔(𝐱,0) = 0 in 𝛺(𝐱,0) (5) where 𝑘𝑑and 𝑘𝑔represent death and growth constants, and the superscripts mean 𝑑(death) and 𝑔(growth), respectively. Finally, the total cell concentration can then be computed as, 𝑛(𝐱, 𝑡) = 𝑛0+𝑛𝑔(𝐱, 𝑡) + 𝑛𝑑(𝐱, 𝑡)(6) Otherwise, Ambrosi and Mollica [23], in place of cell proliferation and death, stated that the increase in the size of the spheroid is due to an increase in the size of the particles that initially make it up. Bull et al. [27], assumed that cell proliferation and death depend on cell cycle times in their agent-based modeling. 3.3. Cell mechanics The spheroid is assumed to be subjected to the contractile activity of cells (see Fig. 1C). This phenomenon is modeled by distinguishing between passive and active components of the cytoskeleton. In particular, contractile forces are implemented as an internal active deformation, i.e. contractility, similar to other existing models of (continuum) cytoskeletal behavior [35,36,57]. For the sake of simplicity of our numerical framework, we will use an updated Lagrangian configuration to solve the equations above. Therefore, a small strains formulation (linearized for time 𝑡around the current configuration 𝛺(𝐱, 𝑡)) is applied to describe the cell mechanical interactions in the spheroid for each small time step. Moreover, we assume a linearly elastic and isotropic constitutive relation between stress and strain around the current configuration. Then, for a material point 𝐱which represents a certain cell concentration, the overall deformation 𝜺𝑇is the contribution of the elastic deformation (𝜀𝑒𝑙) of the passive component of the cell’s cytoskeleton, plus the deformation due to the internal cell contractility 𝜺𝑐𝑒𝑙𝑙. The interpretation of this strain measure is the logarithmic or Hencky strain, so the considered constitutive model is assumed to be dependent on this metric. This is a natural way of describing constitutive stress–strain relations that has been used in a Table 4 Values of the calibrated model parameters and bibliographic range. Parameter Value (U-251 cells) Value (U-87 cells) Bibliographic data Bibliographic data Units 𝑘𝑐0.0011 0.0016 [3.96 ⋅10−8; 0.3] [39,48]mm2 min 𝑟2.9777 ⋅10−10 4.9182 ⋅10−8 [2.4⋅10−15; 0.045] [49,50]mol min⋅cell 𝑘𝑑5⋅10−6 5⋅10−4 [2.73 ⋅10−5; 3.53 ⋅10−4][39,51]mm3 mol⋅min 𝑘𝑔5⋅10−6 5⋅10−4 [8.33 ⋅10−6; 1.034 ⋅10−3][52,53]mm3 mol⋅min 𝑘𝑒𝑔 9.9256 ⋅10−7 1.1476 ⋅10−6 – – mm3 cell 𝑘𝑒𝑑 9.9256 ⋅10−7 8.1969 ⋅10−10 – – mm3 cell 𝑛∗3.0225 ⋅1056.7098 ⋅104– – cell mm3 𝑐𝑛0.5 0.5 – – mol mm3 𝑐𝑠𝑎𝑡 0.9 0.9 – – mol mm3 𝛼4⋅10−8 4⋅10−8 – – mm3 cell⋅min Computers in Biology and Medicine 159 (2023) 106897 6 A. Carrasco-Mantis et al. Fig. 4. Parametric analysis of the distribution of the concentration of living cells along the dimensionless radius of the spheroid, at time steps 0, 1, 3, 5, 7, 9, 11, 13 and 15, for cases 2r, 0.5r, 2k, 0.5k, 2ke, 0.5ke, 2𝑛∗and 0.5𝑛∗. number of hyperelastic models for elastomers and living tissues [58– 60]. Hence, 𝜺𝑇=𝜺𝑒𝑙 +𝜺𝑐𝑒𝑙𝑙 in 𝛺(𝐱, 𝑡)(7) The cell contractility 𝜺𝑐𝑒𝑙𝑙 is assumed to be isotropically dependent on the activity of cells, and hence on cell concentration 𝑛, as follows: 𝜺𝑐𝑒𝑙𝑙 =𝑘𝜀(𝑛)⋅𝑛⋅𝐈in 𝛺(𝐱, 𝑡)(8) where 𝑘𝜀is the contractility constant, which depends on the available room for cells to exert contractile forces. Fig. 3 shows a sketch which explains this hypothesis: given an initial state of contraction 𝑛0[30,37, 38] (red circumference with green cells), two different fates can occur, depending on whether the net balance is growth (G) or death (D). In the case of growth, if there is enough room to accommodate more cells, compression (G.I) occurs because the proliferating cells (green plus purple cells) exert contractile forces. When the reference value is overpassed (value of 𝑛∗is exceeded), there is no more room for cells and the reference volume expands (G.II and G.III). On the other hand, if the fate is death, cells die by necrosis, and then release their contractile forces producing the expansion of the spheroid (D.I). The function which models the cell’s contractility/expansion, according to the hypothesis above, is shown in Fig. 3. Therefore, the contractility constant 𝑘𝜀is defined as, 𝑘𝜀={𝑘𝜀,𝑔 if 𝑛≥𝑛∗ −𝑘𝜀,𝑑 if 𝑛<𝑛∗(9) Several works also consider the shift from compression to expansion. On the one hand, Kopanska et al. [61] studied, from an in vitro point of view, the contractile forces that occur in a spheroid as a consequence of its interaction with the extracellular matrix, concluding that during growth, it undergoes both contraction and expansion as a consequence of its interplay with collagen. On the other hand, Giverso and Preziosi [62] studied the spheroid as a porous material composed of cells and water and established that when it is subjected to radial compression on its surface, it initially undergoes a compressive state with subsequent expansion due to proliferation predominates. In another study, Chen et al. [22], modeled the mechanical interaction between the extracellular matrix and the tumor spheroid, rather than the mechanics among the cells of the spheroid and the tumor spheroid. Mascheroni et al. [19] did not consider a cellular concentration that Computers in Biology and Medicine 159 (2023) 106897 7 A. Carrasco-Mantis et al. Fig. 5. Parametric analysis of the distribution of the concentration of dead cells along the dimensionless radius of the spheroid, at time steps 0, 1, 3, 5, 7, 9, 11, 13 and 15, for cases 2r, 0.5r, 2k, 0.5k, 2ke, 0.5ke, 2𝑛∗and 0.5𝑛∗. marks the transition between zones of spheroidal contraction and expansion. In our model the shift from compression to expansion during growth is aligned with the experimental observed behavior, and is regulated by the model parameter 𝑛∗. The remaining equilibrium, compatibility and constitutive equations (linearized for time 𝑡around the current configuration 𝛺(𝐱, 𝑡)) read as, ∇⋅𝝈=𝟎 𝜺𝑇=1 2(∇𝐮+ ∇𝑇𝐮)in 𝛺(𝐱, 𝑡) 𝝈=𝐂𝜺𝑒𝑙 𝝈(𝐱,0) = 𝝈0in 𝛺(𝐱,0) 𝐭𝑛(𝐱, 𝑡) =  𝐭𝑛in 𝛤𝑡(𝐱, 𝑡) 𝐮(𝐱, 𝑡) =  𝐮in 𝛤𝑢(𝐱, 𝑡) (10) with 𝛤𝑏=𝛤𝑢∪𝛤𝑡,𝛤𝑢∩𝛤𝑡= ∅, and 𝐂being the fourth order elasticity tensor which characterizes the mechanical behavior of the spheroid. 𝝈0 is a certain initial pre-stress in the spheroid, and  𝐭𝑛and  𝐮the prescribed tractions and displacements on the spheroid boundary. These boundary conditions are included in the formulation for completeness purposes, although the analyzed spheroids were considered as free of initial pre-stress, forces and of prescribed displacements along the boundary. 4. Dimensionless formulation Space and time are now non-dimensionalized as follows:  𝒙=𝐱 𝐷0 (11)  𝑡=𝑡 𝑇0 (12) with 𝐷0the diameter of the spheroid, and 𝑇0a characteristic time of analysis. On the other hand, dimensionless nutrient concentration is defined with reference to the initial concentration of nutrients in the suspension experiment: 𝑐( 𝒙, 𝑡) = 𝑐( 𝒙, 𝑡) 𝑐𝛤 (13) while the cell concentration is non-dimensionalized versus the initial cell concentration of the spheroid, 𝑛( 𝒙, 𝑡) = 𝑛( 𝒙, 𝑡) 𝑛0 (14) Computers in Biology and Medicine 159 (2023) 106897 8 A. Carrasco-Mantis et al. Fig. 6. Parametric analysis of the distribution of the concentration of nutrients along the dimensionless radius of the spheroid, at time steps 0, 1, 3, 5, 7, 9, 11, 13 and 15, for cases 2r, 0.5r, 2k, 0.5k, 2ke, 0.5ke, 2𝑛∗and 0.5𝑛∗. We assume a linearly elastic isotropic homogeneous constitutive model for the spheroid, defined by two parameters: the elastic modulus (𝐸) and the Poisson’s ratio (𝜈). The stresses are then referred to the elastic modulus of the spheroid:  𝝈( 𝒙, 𝑡) = 𝝈( 𝒙, 𝑡) 𝐸(15) 𝜈was set to 0.3 in this study. 4.1. Nutrient diffusion Introducing the dimensionless quantities in Eqs. (11)–(14) in the nutrient diffusion Eq. (1) yields, 𝜕 𝑐 𝜕 𝑡= 𝑘𝑐∇2𝑐 −𝑟 ⋅𝑛 in  𝛺( 𝒙, 𝑡) 𝑐( 𝒙,0) = 𝑐0in  𝛺( 𝒙,0) 𝑐( 𝒙, 𝑡) =  𝑓( 𝒙, 𝑡)in  𝛤𝑏( 𝒙, 𝑡) (16) with: 𝜕 𝑓 𝜕𝑡 = −𝛼 ⋅𝑛 in  𝛤𝑏( 𝒙, 𝑡)  𝑓( 𝒙,0) = 1 in  𝛤𝑏( 𝒙,0) (17) where the following dimensionless parameters in Eq. (16) are defined: 𝑐0=𝑐0 𝑐𝛤  𝑘𝑐=𝑘𝑐𝑇0 𝐷2 0 𝑟 =𝑟𝑛0𝑇0 𝑐𝛤 𝛼 =𝛼𝑛0𝑇0 𝑐𝛤 (18) 4.2. Cell dynamics The dimensionless quantities defined in Eqs. (11)–(14) are used in Eqs. (4) and (5), getting, 𝜕 𝑛𝑑 𝜕 𝑡= − 𝑘𝑑⋅(𝑐𝑛−𝑐)⋅𝑛 if 𝑐 ≤𝑐𝑛in  𝛺( 𝒙, 𝑡) 𝑛𝑑( 𝒙,0) = 0 in  𝛺( 𝒙,0) (19) Computers in Biology and Medicine 159 (2023) 106897 9 A. Carrasco-Mantis et al. Fig. 7. Parametric analysis of the distribution of pressure along the dimensionless radius of the spheroid, at time steps 0, 1, 3, 5, 7, 9, 11, 13 and 15, for cases 2r, 0.5r, 2k, 0.5k, 2ke, 0.5ke, 2𝑛∗and 0.5𝑛∗. 𝜕 𝑛𝑔 𝜕 𝑡= 𝑘𝑔⋅(𝑐 −𝑐𝑛)⋅𝑛 if 𝑐 > 𝑐𝑛in  𝛺( 𝒙, 𝑡) 𝜕 𝑛𝑔 𝜕 𝑡= 𝑘𝑔⋅(𝑐𝑠𝑎𝑡 −𝑐𝑛)⋅𝑛 if 𝑐 > 𝑐𝑠𝑎𝑡 in  𝛺( 𝒙, 𝑡) 𝑛𝑔( 𝒙,0) = 0 in  𝛺( 𝒙,0) (20) 𝑛( 𝒙, 𝑡) = 1 + 𝑛𝑔( 𝒙, 𝑡) + 𝑛𝑑( 𝒙, 𝑡)(21) appearing the following dimensionless parameters in Eqs. (19) and (20):  𝑘𝑑=𝑘𝑑𝑇0 𝑛0𝑐𝛤  𝑘𝑔=𝑘𝑔𝑇0 𝑛0𝑐𝛤 𝑐𝑛=𝑐𝑛 𝑐𝛤 𝑐𝑠𝑎𝑡 =𝑐𝑠𝑎𝑡 𝑐𝛤 (22) Now, taking the time derivative of Eq. (21) yields, 𝜕 𝑛 𝜕 𝑡=𝜕 𝑛𝑔 𝜕 𝑡+𝜕 𝑛𝑑 𝜕 𝑡(23) Considering Eqs. (19) and (20), Eq. (23) may be rewritten as: 𝜕 𝑛 𝜕 𝑡= − 𝑘𝑑⋅(𝑐𝑛−𝑐)⋅𝑛 ⋅𝛿1+ 𝑘𝑔⋅(𝑐 −𝑐𝑛)⋅𝑛 ⋅𝛿2 + 𝑘𝑔⋅(𝑐𝑠𝑎𝑡 −𝑐𝑛)⋅𝑛 ⋅𝛿3in  𝛺( 𝒙, 𝑡) 𝑛( 𝒙,0) = 1 in  𝛺( 𝒙,0) (24) with, ⎧ ⎪ ⎨ ⎪ ⎩ 𝛿1= 1 if 𝑐 ≤𝑐𝑛and 𝛿1= 0 elsewhere 𝛿2= 1 if 𝑐𝑛< 𝑐 ≤𝑐𝑠𝑎𝑡 and 𝛿2= 0 elsewhere 𝛿3= 1 if 𝑐 > 𝑐𝑠𝑎𝑡 and 𝛿3= 0 elsewhere (25) Eq. (24) will be used to compute the cell evolution in the spheroid. Dead and alive cells can be computed afterwards from Eqs. (19) and (20). Computers in Biology and Medicine 159 (2023) 106897 16 A. Carrasco-Mantis et al. References [1] Ana S. Nunes, Andreia S. Barros, Elisabete C. Costa, André F. Moreira, Ilídio J. Correia, 3D tumor spheroids as in vitro models to mimic in vivo human solid tumors resistance to therapeutic drugs, Biotechnol. Bioeng. 116 (1) (2019) 206–226. [2] Shawn P. Carey, Alina Starchenko, Alexandra L. McGregor, Cynthia A. Reinhart- King, Leading malignant cells initiate collective epithelial cell invasion in a three-dimensional heterotypic tumor spheroid model, Clin. Exp. Metastasis 30 (5) (2013) 615–630. [3] Aleksandra Karolak, Dmitry A. Markov, Lisa J. McCawley, Katarzyna A. Rejniak, Towards personalized computational oncology: from spatial models of tumour spheroids, to organoids, to tissues, J. R. Soc. Interface 15 (138) (2018) 20170703. [4] Lucia Le Roux, Dawid Schellingerhout, Andrei Volgin, David Maxwell, Katashi Ishihara, Juri Gelovani, Optimizing imaging of 3D multicellular tumor spheroids with fluorescent reporter proteins using confocal microscopy, Microsc. Microanal. 14 (S2) (2008) 734–735. [5] Louis-Bastien Weiswald, Dominique Bellet, Virginie Dangles-Marie, Spherical cancer models in tumor biology, Neoplasia 17 (1) (2015) 1–15. [6] Ibiayi Dagogo-Jack, Alice T. Shaw, Tumour heterogeneity and resistance to cancer therapies, Nat. Rev. Clin. Oncol. 15 (2) (2018) 81. [7] Robert M. Sutherland, Cell and environment interactions in tumor microregions: the multicell spheroid model, Science 240 (4849) (1988) 177–184. [8] Juergen Friedrich, Claudia Seidel, Reinhard Ebner, Leoni A. Kunz-Schughart, Spheroid-based drug screen: considerations and practical approach, Nat. Protoc. 4 (3) (2009) 309. [9] Sritama Nath, Gayathri R. Devi, Three-dimensional culture systems in cancer research: Focus on tumor spheroid model, Pharmacol. Ther. 163 (2016) 94–108. [10] Christoph Mark, Thomas J. Grundy, Pamela L. Strissel, David Böhringer, Nadine Grummel, Richard Gerum, Julian Steinwachs, Carolin C. Hack, Matthias W. Beckmann, Markus Eckstein, et al., Collective forces of tumor spheroids in three-dimensional biopolymer networks, Elife 9 (2020) e51912. [11] Jose M. Ayuso, Haneen A. Basheer, Rosa Monge, Pablo Sánchez-Álvarez, Manuel Doblaré, Steven D. Shnyder, Victoria Vinader, Kamyar Afarinkia, Luis J. Fernández, Ignacio Ochoa, Study of the chemotactic response of multicellular spheroids in a microfluidic device, PLoS One 10 (10) (2015) e0139515. [12] Ludivine Guillaume, Lise Rigal, Jérôme Fehrenbach, Childérick Severac, Bernard Ducommun, Valérie Lobjois, Characterization of the physical properties of tumorderived spheroids reveals critical insights for pre-clinical studies, Sci. Rep. 9 (1) (2019) 1–9. [13] Jiandong Wu, Xun Wu, Francis Lin, Recent developments in microfluidics-based chemotaxis studies, Lab Chip 13 (13) (2013) 2484–2499. [14] Shing-Yi Cheng, Steven Heilman, Max Wasserman, Shivaun Archer, Michael L. Shuler, Mingming Wu, A hydrogel-based microfluidic device for the studies of directed cell migration, Lab Chip 7 (6) (2007) 763–769. [15] Jose M. Ayuso, Rosa Monge, Guillermo A. Llamazares, Marco Moreno, Maria Agirregabiria, Javier Berganzo, Manuel Doblare, Iñaki Ochoa, Luis J. Fernandez, SU-8 based microdevices to study self-induced chemotaxis in 3D microenvironments, Front. Mater. 2 (2015) 37. [16] Jose M. Ayuso, María Virumbrales-Muñoz, Alodia Lacueva, Pilar M. Lanuza, Elisa Checa-Chavarria, Pablo Botella, Eduardo Fernández, Manuel Doblare, Simon J. Allison, Roger M. Phillips, et al., Development and characterization of a microfluidic model of the tumour microenvironment, Sci. Rep. 6 (1) (2016) 1–16. [17] Aliya Fatehullah, Si Hui Tan, Nick Barker, Organoids as an in vitro model of human development and disease, Nature Cell Biol. 18 (3) (2016) 246–254. [18] Shi-Lei Xue, Bo Li, Xi-Qiao Feng, Huajian Gao, Biochemomechanical poroelastic theory of avascular tumor growth, J. Mech. Phys. Solids 94 (2016) 409–432. [19] Pietro Mascheroni, Cinzia Stigliano, Melania Carfagna, Daniela P. Boso, Luigi Preziosi, Paolo Decuzzi, Bernhard A. Schrefler, Predicting the growth of glioblastoma multiforme spheroids using a multiphase porous media model, Biomech. Model. Mechanobiol. 15 (5) (2016) 1215–1228. [20] H.P. Greenspan, On the growth and stability of cell cultures and solid tumors, J. Theoret. Biol. 56 (1) (1976) 229–242. [21] H.M. Byrne, Mark A.J. Chaplain, Free boundary value problems associated with the growth and development of multicellular spheroids, European J. Appl. Math. 8 (6) (1997) 639–658. [22] C.Y. Chen, H.M. Byrne, J.R. King, The influence of growth-induced stress from the surrounding medium on the development of multicell spheroids, J. Math. Biol. 43 (3) (2001) 191–220. [23] D. Ambrosi, F. Mollica, On the mechanics of a growing tumor, Internat. J. Engrg. Sci. 40 (12) (2002) 1297–1316. [24] D. Ambrosi, F. Mollica, The role of stress in the growth of a multicell spheroid, J. Math. Biol. 48 (5) (2004) 477–499. [25] D. Ambrosi, Simone Pezzuto, Davide Riccobelli, T. Stylianopoulos, P. Ciarletta, Solid tumors are poroelastic solids with a chemo-mechanical feedback on growth, J. Elasticity 129 (1) (2017) 107–124. [26] Massimiliano Fraldi, Angelo R. Carotenuto, Cells competition in tumor growth poroelasticity, J. Mech. Phys. Solids 112 (2018) 345–367. [27] Joshua A. Bull, Franziska Mech, Tom Quaiser, Sarah L. Waters, Helen M. Byrne, Mathematical modelling reveals cellular dynamics within tumour spheroids, PLoS Comput. Biol. 16 (8) (2020) e1007961. [28] M. Amereh, R. Edwards, M. Akbari, B. Nadler, In-silico modeling of tumor spheroid formation and growth, Micromachines 12 (2021) 749. [29] Holger Perfahl, Barry D. Hughes, Tomás Alarcón, Philip K. Maini, Mark C. Lloyd, Matthias Reuss, Helen M. Byrne, 3D hybrid modelling of vascular network formation, J. Theoret. Biol. 414 (2017) 254–268. [30] Jacobo Ayensa-Jiménez, Marina Pérez-Aliacar, Teodora Randelovic, Sara Oliván, Luis Fernández, José Antonio Sanz-Herrera, Ignacio Ochoa, Mohamed H. Doweidar, Manuel Doblaré, Mathematical formulation and parametric analysis of in vitro cell models in microfluidic devices: application to different stages of glioblastoma evolution, Sci. Rep. 10 (1) (2020) 1–21. [31] Anna Kane Laird, Dynamics of tumour growth, Br. J. Cancer 18 (3) (1964) 490. [32] Graeme J. Pettet, C.P. Please, M.J. Tindall, D.L.S. McElwain, The migration of cells in multicell tumor spheroids, Bull. Math. Biol. 63 (2) (2001) 231–257. [33] Xin Lu, Yibin Kang, Hypoxia and hypoxia-inducible factors: master regulators of metastasis, Clin. Cancer Res. 16 (24) (2010) 5928–5935. [34] A.S.G. Curtis, G.M. Seehar, The control of cell division by tension or diffusion, Nature 274 (5666) (1978) 52–53. [35] P. Moreo, J.M. Garcia-Aznar, M. Doblare, Modeling mechanosensing and its effect on the migration and proliferation of adherent cells, Acta Biomater. 4 (2008) 613–621. [36] J.A. Sanz-Herrera, P. Moreo, J.M. Garcia-Aznar, M. Doblare, On the effect of substrate curvature on cell mechanics, Biomaterials 30 (34) (2009) 6674–6686. [37] Sui Huang, Donald E. Ingber, Cell tension, matrix mechanics, and cancer development, Cancer Cell 8 (3) (2005) 175–176. [38] Anika Nagelkerke, Johan Bussink, Alan E. Rowan, Paul N. Span, The mechanical microenvironment in cancer: How physics affects tumours, in: Seminars in Cancer Biology, 35, Elsevier, 2015, pp. 62–70. [39] Alicia Martínez-González, Gabriel F. Calvo, Luis A. Pérez Romasanta, Víctor M. Pérez-García, Hypoxic cell waves around necrotic cores in glioblastoma: a biomathematical model and its therapeutic implications, Bull. Math. Biol. 74 (12) (2012) 2875–2896. [40] Daniel J. Brat, Glioblastoma: biology, genetics, and behavior, in: American Society of Clinical Oncology Educational Book, Vol. 32, No. 1, American Society of Clinical Oncology, 2012, pp. 102–107. [41] David Robert Grimes, Catherine Kelly, Katarzyna Bloch, Mike Partridge, A method for estimating the oxygen consumption rate in multicellular tumour spheroids, J. R. Soc. Interface 11 (92) (2014) 20131124. [42] Hongmei Yu, Janna Kay Mouw, Valerie M. Weaver, Forcing form and function: biomechanical regulation of tumor evolution, Trends Cell Biol. 21 (1) (2011) 47–56. [43] David Lacalle, Héctor Alfonso Castro-Abril, Teodora Randelovic, César Domínguez, Jónathan Heras, Eloy Mata, Gadea Mata, Yolanda Méndez, Vico Pascual, Ignacio Ochoa, SpheroidJ: An open-source set of tools for spheroid segmentation, Comput. Methods Programs Biomed. 200 (2021) 105837. [44] Priya R. Baraniak, Todd C. McDevitt, Scaffold-free culture of mesenchymal stem cell spheroids in suspension preserves multilineage potential, Cell Tissue Res. 347 (2012) 701–711. [45] Ian F. Tannock, Oxygen diffusion and the distribution of cellular radiosensitivity in tumours, Br. J. Radiol. 45 (535) (1972) 515–524. [46] Brian W. Pogue, Julia A. O’Hara, Carmen M. Wilmot, Keith D. Paulsen, Harold M. Swartz, Estimation of oxygen distribution in RIF-1 tumors by diffusion modelbased interpretation of pimonidazole hypoxia and Eppendorf measurements, Radiat. Res. 155 (1) (2001) 15–25. [47] Peggy L. Olive, Charlene Vikse, Martin J. Trotter, Measurement of oxygen diffusion distance in tumor cubes using a fluorescent hypoxia probe, Int. J. Radiat. Oncol. Biol. Phys. 22 (3) (1992) 397–402. [48] Christina H. Wang, Jason K. Rockhill, Maciej Mrugala, Danielle L. Peacock, Albert Lai, Katy Jusenius, Joanna M. Wardlaw, Timothy Cloughesy, Alexander M. Spence, Russ Rockne, et al., Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model, Cancer Res. 69 (23) (2009) 9133–9140. [49] Abramo Agosti, Clara Cattaneo, Chiara Giverso, Davide Ambrosi, Pasquale Ciarletta, A computational framework for the personalized clinical treatment of glioblastoma multiforme, ZAMM Z. Angew. Math. Mech. 98 (12) (2018) 2307–2327. [50] Philip Gerlee, Alexander R.A. Anderson, An evolutionary hybrid cellular automaton model of solid tumour growth, J. Theoret. Biol. 246 (4) (2007) 583–603. [51] Abramo Agosti, Chiara Giverso, Elena Faggiano, Aymeric Stamm, Pasquale Ciarletta, A personalized mathematical tool for neuro-oncology: A clinical case study, Int. J. Non-Linear Mech. 107 (2018) 170–181. [52] Corinne E. Griguer, Claudia R. Oliva, Eric Gobin, Pascale Marcorelles, Dale J. Benos, Jack R. Lancaster Jr., G. Yancey Gillespie, CD133 is a marker of bioenergetic stress in human glioma, PLoS One 3 (11) (2008) e3655. [53] Alicia Martínez-González, Mario Durán-Prado, Gabriel F. Calvo, Francisco J. Alcaín, Luis A. Pérez-Romasanta, Víctor M. Pérez-García, Combined therapies of antithrombotics and antioxidants delay in silico brain tumour progression, Math. Med. Biol. 32 (3) (2015) 239–262. Computers in Biology and Medicine 159 (2023) 106897 17 A. Carrasco-Mantis et al. [54] H.M. Byrne, T. Alarcon, M.R. Owen, S.D. Webb, P.K. Maini, Modelling aspects of cancer dynamics: a review, Phil. Trans. R. Soc. A 364 (1843) (2006) 1563–1578. [55] Jane S. Tjia, Prabhas V. Moghe, Regulation of cell motility on polymer substrates via" dynamic," cell internalizable, ligand microinterfaces, Tissue Eng. 8 (2) (2002) 247–261. [56] Brian Stramer, Roberto Mayor, Mechanisms and in vivo functions of contact inhibition of locomotion, Nature Rev. Mol. Cell. Biol. 18 (1) (2017) 43–55. [57] M.R.K. Mofrad, R.D. Kamm, Continuum Elastic or Viscoelastic Models for the Cell, Cambridge University Press, 2006, [58] T. Sussman, K.J. Bathe, A model of incompressible isotropic hyperelastic material behavior using spline interpolations of tension–compression test data, Commun. Numer. Meth. Eng. 25 (2009) 53–63. [59] M. Latorre, F.J. Montans, Extension of the Sussman–Bathe spline-based hyperelastic model to incompressible transversely isotropic materials, Comput. Struct. 122 (2013) 13–26. [60] M. Latorre, F.J. Montans, On the interpretation of the logarithmic strain tensor in an arbitrary system of representation, Int. J. Solids Struct. 51 (2014) 1507–1515. [61] Katarzyna S. Kopanska, Yara Alcheikh, Ralitza Staneva, Danijela Vignjevic, Timo Betz, Tensile forces originating from cancer spheroids facilitate tumor invasion, PLoS One 11 (6) (2016) e0156442. [62] Chiara Giverso, Luigi Preziosi, Influence of the mechanical properties of the necrotic core on the growth and remodelling of tumour spheroids, Int. J. Non-Linear Mech. 108 (2019) 20–32. [63] O.C. Zienkiewicz, Robert L. (Robert Leroy) Taylor, The Finite Element Method, Butterworth–Heinemann, 2000. [64] J.A. Sanz-Herrera, J. Barrasa-Fano, M. Condor, H. Van Oosterwyck, Inverse method based on 3D nonlinear physically constrained minimisation in the framework of traction force microscopy, Soft Matter 17 (2021) 10210–10222. [65] J. Barrasa-Fano, A. Shapeti, J. De Jong, A. Ranga, H. Sanz-Herrera, Advanced in silico validation framework for three-dimensional traction force microscopy and application to an in vitro model of sprouting angiogenesis, Acta Biomater. 126 (2021) 326–338. [66] Richard M. Enmon Jr., Kim C. O’Connor, Daniel J. Lacks, Daniel K. Schwartz, Robert S. Dotson, Dynamics of spheroid self-assembly in liquid-overlay culture of DU 145 human prostate cancer cells, Biotechnol. Bioeng. 72 (6) (2001) 579–591. [67] Meitham Amereh, Roderick Edwards, Mohsen Akbari, Ben Nadler, In-silico modeling of tumor spheroid formation and growth, Micromachines 12 (7) (2021) 749. [68] I. Smyrek, B. Mathew, S.C. Fischer, S.M. Lissek, S. Becker, E.H.K. Stelzer, E- cadherin, actin, microtubules and FAK dominate different spheroid formation phases and important elements of tissue integrity, Biol. Open 8 (1) (2019) bio037051. [69] Stephanie Eugenie Brigitte McArdle, Kinana Habra, Joshua R.D. Pearson, Swift formation of optimal single spheroids towards in-vitro 3-dimensional tumour models, 2021, bioRxiv. [70] X. Cui, Y. Hartanto, H. Zhang, Advances in multicellular spheroids formation, J. R. Soc. Interface 14 (127) (2017) 20160877. [71] Ruei-Zeng Lin, Li-Fang Chou, Chi-Chen Michael Chien, Hwan-You Chang, Dynamic analysis of hepatoma spheroid formation: roles of E-cadherin and 𝛽1-integrin, Cell Tissue Res. 324 (3) (2006) 411–422. [72] Jiangping Xu, Guillermo Vilanova, Hector Gomez, A mathematical model coupling tumor growth and angiogenesis, PLoS One 11 (2) (2016) e0149422. [73] Abdullah Chandra Sekhar Talari, Ahtasham Raza, Shazza Rehman, Ihtesham U. Rehman, Analyzing normal proliferating, hypoxic and necrotic regions of T-47D human breast cancer spheroids using Raman spectroscopy, Appl. Spectrosc. Rev. 52 (10) (2017) 909–924. [74] Hezhen Li, Bingxi Lei, Wei Xiang, Hai Wang, Wenfeng Feng, Yawei Liu, Songtao Qi, Differences in protein expression between the U251 and U87 cell lines, Turk Neurosurg 27 (6) (2017) 894–903.