Examining avascular tumour growth dynamics: A variable-order non-local modelling perspective

Abstract

This study investigates the growth of an avascular tumour described through the interchange of mass among its constituents and the production of inelastic distortions induced by growth itself. A key contribution of this research examines the role of non-local diffusion arising from the complex and heterogeneous tumour micro-environment. In our context, the non-local diffusion is enhanced by a variable-order fractional operator that incorporates crucial information about regions of limited nutrient availability within the tissue. Our research also focuses on the identification of an evolution law for the growth-induced inelastic distortions recast through the identification of non-conventional forces dual to suitable kinematic descriptors associated with the growth tensor. The establishment of such evolution law stems from examining the dissipation inequality and subsequently determining a posteriori connections between the inelastic distortions and the source/sink terms in the mass balance laws. To gain insights into the dynamics of tumour growth and its response to the proposed modelling framework, we first study how the variables governing the tissue evolution are affected by the introduction of the new growth law. Second, we investigate how regions of limited diffusion within the tumour, encoded into a fractional operator of variable-order, influence its growth.

1. Introduction

Cancer is a complex disease and a prominent cause of mortality worldwide [1] with one common aspect being its diversity [2]. The inherent complexity of cancer poses significant challenges for studying and comprehending the disease and complicates the discovery of treatments that can be effective for all patients [2, 3]. Owing to its intricate nature and the diverse range of diseases it encompasses, cancer lacks a universal description, making its aetiology a subject of ongoing debate [3–5]. Nevertheless, there is a general consensus that both inherited and external factors contribute to the development of mutations in gene expression transforming normal cells into tumour cells [1, 4–7]. A multistage process characterised by an imbalance between cell proliferation and cell death leads to the emergence of a population of cells capable of invasive growth and metastasis [1, 6, 8]. Tumour progression can be categorised into three distinct stages: avascular, vascular, and metastatic [7, 9, 10].

Avascular tumours, also known as in situ cancers, lack their own blood vessels and rely on diffusion for oxygen and nutrients supply [4, 10]. This limited blood supply allows the tumour to grow slowly and remain dormant for extended periods. However, it can transition into a dangerous vascular tumour through angiogenesis [7]. In laboratory settings, the initial avascular phase of tumour development can be studied using three-dimensional multicellular spheroids composed of cancer cells [7, 11–13]. In vitro experiments using tumour spheroids have shown that their growth ceases once they reach a specific size, because of the balance between cell proliferation and necrosis [12]. Proliferating cells are primarily found in the outer layers, necrotic cells in the centre, and quiescent cells in an intermediate environment, with some capable of being recruited to the proliferating cell layer [7, 12, 14].

Considerable efforts have been devoted to investigating tumour progression through extensive research (see, for instance, [7, 9, 13–26]). These models provide a means to explore and predict tumour behaviour that may not be readily observable in experimental studies, potentially leading to novel therapies or improvements in existing treatments [23]. However, owing to the complexity of tumour behaviour and limited knowledge of growth processes [2–5], the development of realistic models remains a significant challenge [27,28]. Mathematical models describing tumour growth can be classified into discrete and continuum models. A comprehensive review of these models can be found in Roose et al. [7]. In this study, we focus on utilising a continuum approach to model the growth of avascular tumours. Although modelling avascular solid tumour growth might not be as essential for cancer therapy as modelling vascular and angiogenic tumour growth, understanding the evolution of avascular tumours can provide insights into more complex models. Avascular tumour growth is simpler to model mathematically, has abundant experimental data, and encompasses concepts applicable to general vascular tumour growth models [7]. Continuum models account for the interaction between cell density and chemical agents, such as nutrients, which influence tumour cell cycle events and provide insights into tumour growth dynamics [7]. These models are often based on the Theory of Mixtures [29–35] where tumour growth is accompanied by inelastic distortions and can be associated with chemo-mechanical interactions [17, 36, 37]. In this respect, the Bilby–Kroner–Lee (BKL) decomposition of the deformation gradient tensor [37–40] has been commonly used to account for alterations in the internal structure of the tissue caused by changes in mass.

Nutrients play a significant role in tumour progression because cancer cells rely on them for proliferation [7, 12, 14, 41, 42]. Fick’s law–based diffusion models are commonly employed to describe the transport of nutrients. However, recent studies have identified anomalous diffusion patterns, also known as non-Fickian or non-local, in biological tissues [43–46]. Non-local diffusion, which considers interactions between distant particles, has been attributed to the multiscale and heterogeneous nature of the medium in which it is taking place [44, 46, 47]. For instance, as discussed in Höfling and Franosch [45], macromolecular crowding leads to densely packed and heterogeneous structures, causing subdiffusive behaviour in the mean-square displacement and spatially non-Gaussian particle movement, thereby challenging the standard Brownian motion assumption.

Fractional calculus [48–50] offers a mathematical framework for investigating non-Fickian transport processes and anomalous diffusion, and it has shown promise in understanding cancer-related phenomena such as chemotherapy, radiotherapy, and tumour growth [26, 51–54]. The utilisation of integro-differential operators, characterised by a constant order $\alpha \in \mathbb{R}$ , has proven to be highly effective in explaining diverse non-local phenomena [47, 55–57]. However, there is a rising interest in fractional operators that incorporate a variable-order (depending on time, space or specific variables) for representing evolutionary phenomena while keeping the governing equations unchanged [58]. Notably, variable-order integrodifferential operators are increasingly recognised as valuable tools in modelling complex real-world problems [58–61]. These include the study of the mechanics of materials [62–64], the dynamics of biological tissues [25, 65], transport processes and diffusion [66–68], among others. These examples demonstrate the potential of variable-order fractional operators in advancing our comprehension of various fields, enabling more precise and comprehensive modelling of intricate systems that exhibit non-local effects.

In the present work, our objective is to examine the influence of non-local diffusion mechanisms on the growth of an avascular tumour. To achieve this, we describe the tumour tissue as a biphasic medium comprising a fluid and a solid phase, and we associate its growth with the gain or loss of mass of the solid phase at the expense or advantage of the fluid one [9, 24, 26]. In this respect, it is important to note that our current model focuses on the growth of a tumour tissue in the absence of blood vessels, considering only proliferating and necrotic cells in the solid phase and neglecting quiescent cells in our analysis. Furthermore, to describe the diffusion of chemical agents, we slightly modify the constitutive law proposed in Ramírez-Torres et al. [26] by introducing a spatial dependency in the order of the integro-differential operator. This represents an additional feature with respect to the definition given in Ramírez-Torres et al. [26] and will allow us to capture potential spatial phenomena that may arise due to the heterogeneous and complex environment in which diffusion occurs. Our goal is to include relevant information about the complex tumour micro-environment, which can lead to regions with limited difussion [69]. We remark that, although we refer to a variable-order fractional operator in space, our system evolves over time as the variable-order parameter is influenced by the motion of the body, as well as the inelastic distortions it undergoes.

Another extra feature with respect to Ramírez-Torres et al. [26] is that we adopt the methodology proposed in Grillo and Di Stefano [70, 71] to derive a dynamic equation for the evolution of growth-induced inelastic distortions. This methodology relies on the admission of the existence of non-conventional forces [36, 70] whose balance represents the evolution equation for the inelastic distortions. Furthermore, it avoids making initial phenomenological assumptions in the source/sink term corresponding to the mass balance of the solid phase as a whole. Within this setting, the equation for the inelastic distortions is identified through an a posteriori methodology [70, 71], i.e., after the constitutive law governing the internal non-conventional force has been established through the study of the dissipation inequality. In our context, a consequence of this approach involves the necessity of introducing source and sink terms accounting for the influence of configurational mechanical stress in the mass balance equations for the proliferating and necrotic cells. Indeed, once the equation for the growth-induced inelastic distortions has been established, these new terms are attributed to the presence of the Eshelby stress tensor in the constitutive law governing the internal non-conventional force. In particular, they characterise the influence of the inelastic distortions on the dynamics of the solid phase’s constituents. We mention that the introduction of source/sink terms accounting for the impact of inelastic distortions on the growth of a tumour has also been investigated in Di Stefano et al. [24]. Finally, in agreement with previous research [9, 24, 26], we formalise the growth law by expressing the external non-conventional force in such a way it accounts for the progression of nutrients and mechanotransduction.

We conclude this work by examining the consequences of our modelling approach, which aims to understand how the dynamics of an avascular tumour are affected by a non-Fickian diffusion law of variable-order. Our analysis focuses on two main aspects: first, we investigate how the variables governing the tumour’s evolution are influenced by the new growth law, which accounts for configurational mechanical stress, and second, we explore the impact of a specific non-uniform fractional-order parameter encoding information about regions with limited nutrient availability within the tumour tissue. To this end, we investigate a benchmark problem that represents the progression of ductal carcinoma in situ (DCIS). DCIS is increasingly recognised as a target for preventive measures [72] as it is considered a precursor to invasive breast cancer and shares many biological characteristics with invasive diseases. Therefore, the study of DCIS is relevant since, if left untreated, DCIS can potentially develop and infiltrate the breast stromal tissue surrounding the ducts, thereby becoming life-threatening [72, 73].

2. Description of the avascular tumour as a mixture

Here, we briefly introduce some of the main definitions used in the description of a mixture consisting of two phases: a fluid phase, denoted with ${\mathcal{P}}_{\mathrm{f}}$ , and a solid phase ${\mathcal{P}}_{\mathrm{s}}$ [32–35]. Specifically, by considering some of the assumptions made in Mascheroni et al. [9], Di Stefano et al. [24], Ramírez-Torres et al. [26], Mascheroni et al. [74], and Grillo et al. [75] in the case of an avascular tumour, we assume that the solid phase consists of proliferating and necrotic cells, with mass fractions ${\omega }_{\mathrm{ps}}$ and ${\omega }_{\mathrm{ns}}$ , respectively, while the fluid phase is made of water and chemical agents, with mass fractions ${\omega }_{\mathrm{wf}}$ and ${\omega }_{\mathrm{cf}}$ , respectively. In our context, we speak of nutrients instead of chemical agents, and one example is glucose, which serves as an energy source facilitating cell proliferation. Furthermore, we suppose that both the solid and the fluid phases are saturated, i.e., ${\omega }_{\mathrm{ps}}+{\omega }_{\mathrm{ns}}=1$ and ${\omega }_{\mathrm{wf}}+{\omega }_{\mathrm{cf}}=1$ . In addition, the apparent mass densities of the solid and fluid phases are denoted by ${\phi }_{\mathrm{s}}{\varrho }_{\mathrm{s}}$ and ${\phi }_{\mathrm{f}}{\varrho }_{\mathrm{f}}$ , respectively, where ${\varrho }_{\mathrm{s}}$ and ${\varrho }_{\mathrm{f}}$ represent the true mass densities of the solid and fluid phases, and ${\phi }_{\mathrm{s}}$ and ${\phi }_{\mathrm{f}}$ are the solid and the fluid phases volumetric fractions. Here, we also assume to be in the presence of a saturated mixture, so that the condition ${\phi }_{\mathrm{f}}+{\phi }_{\mathrm{s}}=1$ holds true.

2.1. Kinematics

As in previous works [9, 24, 26, 75], we let $\mathcal{S}$ be the three-dimensional Euclidean space and denote with $\mathcal{B}$ and ${\mathcal{B}}_t$ , respectively, the reference and current configurations of the continuum representing the avascular tumour. Following the discussion in Quiligotti et al. [76] and Crevacore et al. [77], ${\mathcal{B}}_t$ corresponds to the region within the Euclidean space $\mathcal{S}$ , where both solid and fluid phases coexist, namely, the solid-fluid mixture. The map $\chi (.,t):\mathcal{B}\to \mathcal{S}$ denotes the motion of the solid phase, and we assume that, for each time $t\in \mathrm{S}$ , the inverse mapping ${\chi }^{-1}(.,t):{\mathcal{B}}_t\to \mathcal{S}$ maps every point of ${\mathcal{B}}_t$ to the reference configuration of the solid phase. Consequently, since ${\chi }^{-1}(.,t)$ is surjective, χ characterises the motion of the entire mixture and $\mathcal{B}$ can be taken as the reference configuration of the mixture. Furthermore, if v _s and v _f are the spatial velocities of the solid and the fluid phase, we can write v _s(x, t) = ${\hat{v}}_{\mathrm{s}}(X,t):=\stackrel{·}{\chi }(X,t)$ and ${\mathit{v}}_{\mathrm{f}}(x,t)={\hat{\mathit{v}}}_{\mathrm{f}}(X,t)$ , where $X\in \mathcal{B}$ and the “hat” notation is used to represent the material form of the spatial field. That is, for a generic spatial field $\varphi$ , $\hat{\varphi }:=\varphi °(\chi ,\mathcal{T})$ , where $\mathcal{T}:\mathcal{B}\times \mathcal{I}\to \mathcal{I}$ is such that $\mathcal{T}(X,t)=t$ , with $\mathcal{I}$ denoting a bounded time interval.

The visible deformation of the mixture is characterised by the deformation gradient tensor F , which is defined as the tangent map of the motion χ and its determinant, denoted by J := det F , is assumed to be positive. Here, we focus on volumetric growth and as in previous works (see, for instance, [37, 39, 78–81]), we describe it as a phenomenon that induces inelastic distortions. In this respect, we introduce the BKL decomposition of the deformation gradient [37–39, 78, 79, 81, 82]:

$$\mathit{F}(X,t)={\mathit{F}}_e(X,t){\mathit{F}}_{\gamma }(X,t),$$

(1)

where the second-order tensor F _e represents the tensor of elastic distortions or elastic accommodation, while ${\mathit{F}}_{\gamma }$ is known as the tensor of inelastic distortions or growth tensor. The latter maps, for each pair $(X,t)\in \mathcal{B}\times \mathcal{I}$ , vectors from the tangent space $T_X\mathcal{B}$ to vectors in ${\mathcal{N}}_t\left(X\right)$ , often referred to as the body’s natural state. More precisely, ${\mathcal{N}}_t\left(X\right)$ represents the state where the tissue is entirely free of stress. If growth is seen as an inelastic process that generates stresses that cannot be eliminated by just removing applied loads, achieving this stress-free state requires an ideal tearing process that effectively disassembles the tissue into a collection of completely relaxed pieces. Experiments conducted on tumours (see, for instance, [13]) support the idea of growth-induced inelastic distortions. Further discussions about the multiplicative decomposition (1) can be found in Di Stefano et al. [24], Goriely [37], Mićnović [38], Sadik and Yavari [82], and Ciancio et al. [83].

We indicate by J_e and $J_{\gamma }$ the determinants of F _e and ${\mathit{F}}_{\gamma }$ , respectively, so that from equation (1), we can write $J=J_{\mathrm{e}}{J}_{\gamma }$ . Furthermore, we notice that the time derivative of $J_{\gamma }$ is:

$${\stackrel{·}{J}}_{\gamma }=J_{\gamma }\mathrm{tr}\left({\mathit{F}}_{\gamma }^{-1}{\stackrel{·}{\mathit{F}}}_{\gamma }\right).$$

(2)

In this work, we build upon the consideration that ${\mathit{F}}_{\gamma }$ is spherical and set ${\mathit{F}}_{\gamma }=\gamma \mathit{I}$ , where γ > 0 is the growth parameter and I denotes the material identity tensor. In particular, this consideration implies that equation (2) can be rewritten in the form:

$$J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }=3{\gamma }^{-1}\stackrel{·}{\gamma }.$$

(3)

3. Balance laws

3.1. Mass balance laws

If the mass of the mixture’s constituents is not conserved, then mass exchange exists among them. In this context, the mass balance laws for the constituents of the mixture can be written as [9, 24, 26, 35, 80, 84]:

$${\partial }_t\left({\phi }_{\mathrm{s}}{\varrho }_{\mathrm{s}}{\omega }_{\mathrm{ps}}\right)+\mathrm{div}\left({\phi }_{\mathrm{s}}{\varrho }_{\mathrm{s}}{\omega }_{\mathrm{ps}}{\mathit{v}}_{\mathrm{s}}\right)=r_{\mathrm{pn}}+r_{\mathrm{fp}}+r_{\mathrm{ph}},$$

(4a)

$${\partial }_t\left({\phi }_{\mathrm{s}}{\varrho }_{\mathrm{s}}\right)+\mathrm{div}\left({\phi }_{\mathrm{s}}{\varrho }_{\mathrm{s}}{\omega }_{\mathrm{ps}}{\mathit{v}}_{\mathrm{s}}\right)=r_{\mathrm{s}},$$

(4b)

$${\partial }_t\left({\phi }_{\mathrm{f}}{\varrho }_{\mathrm{f}}{\omega }_{\mathrm{cf}}\right)+\mathrm{div}({\phi }_{\mathrm{f}}{\varrho }_{\mathrm{f}}{\omega }_{\mathrm{cf}}{\mathit{v}}_{\mathrm{f}}+{\mathit{y}}_{\alpha })=r_{\mathrm{cf}},$$

(4c)

$${\partial }_t\left({\phi }_{\mathrm{f}}{\varrho }_{\mathrm{f}}\right)+\mathrm{div}\left({\phi }_{\mathrm{f}}{\varrho }_{\mathrm{f}}{\omega }_{\mathrm{cf}}{\mathit{v}}_{\mathrm{f}}\right)=r_{\mathrm{f}},$$

(4d)

where equations (4a) and (4b) represent the mass balance for the proliferating cells and the solid phase as a whole, whereas equations (4c) and (4d) characterise the mass balance for the nutrients and the fluid phase as a whole, respectively.

In writing equations (4a)–(4d), we have assumed that the proliferating cells move with the same velocity of the solid phase [9,24,26] and we have denoted by ${\mathit{y}}_{\alpha }$ the mass flux vector of the chemical agents in the fluid phase. We notice that the presence of the subscript α in the notation of the mass flux vector is due to our efforts in deviating from Fick’s law to describe the diffusion of nutrients. To this end, as anticipated in section 1, we benefit from the notion of fractional differentiation [48–50] which considers integro-differential operators of order $\alpha \in \mathbb{R}$ . Although we do not provide at this point a constitutive expression for the mass flux vector ${\mathit{y}}_{\alpha }$ , comprehensive insights into the implications of this selection will be provided in the following sections.

To account for the mass exchange among the constituents of the system under consideration, we have introduced source/sink terms. These are $r_{\mathrm{pn}}$ , which is used to describe the rate at which the proliferating cells die, $r_{\mathrm{fp}}$ takes into consideration the mass originating from the fluid phase that facilitates cell proliferation, and $r_{\mathrm{cf}}$ is the rate at which the proliferating cells consume nutrients. In our formulation, we have also introduced the notation $r_{\mathrm{s}}:=r_{\mathrm{fp}}+r_{\mathrm{nf}}+r_{\mathrm{ph}}+r_{\mathrm{nh}}$ where the rate at which the fluid dissolves necrotic cells is denoted as $r_{\mathrm{nf}}$ and, to ensure mass conservation of the mixture as a whole, we have set $r_{\mathrm{f}}=-r_{\mathrm{s}}$ . We remark that a difference with respect to the model presented in Ramírez-Torres et al. [26] relies on the introduction of the source/sink terms $r_{\mathrm{ph}}$ and $r_{\mathrm{nh}}$ into the balance equations for the proliferating and necrotic cells which, in consonance with Di Stefano et al. [24], characterise potential effects of the properties of ${\mathit{F}}_{\gamma }$ on the growth process. Specifically, as it will be clear in the following sections, $r_{\mathrm{ph}}$ and $r_{\mathrm{nh}}$ enclose the effects of the inelastic distortions induced by growth. Their existence will showcase the contribution of configurational mechanical stress to the dynamics of the tissue and complement those of different nature (e.g., of chemical nature) specified through the external non-conventional force.

As in Giverso et al. [20], Di Stefano et al. [24], Ramírez-Torres et al. [26], and Byrne and Preziosi [85], we consider the true mass densities to be constant and equal for both the solid and fluid phases, namely, ${\varrho }_{\mathrm{s}}(x,t)={\varrho }_{\mathrm{f}}(x,t)\equiv {\varrho }_{\mathrm{s}0}$ . Consequently, the mass balance equations (4a)–(4d), expressed in relation to the reference configuration, are given by:

$${\stackrel{·}{\hat{\omega }}}_{\mathrm{ps}}=J[{\hat{r}}_{\mathrm{pn}}+{\hat{r}}_{\mathrm{fp}}+{\hat{r}}_{\mathrm{ph}}-{\hat{r}}_{\mathrm{s}}{\hat{\omega }}_{\mathrm{ps}}]{\left[J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}\right]}^{-1},$$

(5a)

$${\stackrel{·}{J}}_{\gamma }=J_{\gamma }{R}_{\gamma },$$

(5b)

$${\varrho }_{\mathrm{s}0}[J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}]{\stackrel{·}{\hat{\omega }}}_{\mathrm{cf}}+{\varrho }_{\mathrm{s}0}\mathit{Q}\mathrm{Grad}\left({\hat{\omega }}_{\mathrm{cf}}\right)+\mathrm{Div}\left({\mathit{Y}}_{\alpha }\right)=J[{\hat{r}}_{\mathrm{cf}}+{\hat{r}}_{\mathrm{s}}{\hat{w}}_{\mathrm{cf}}],$$

(5c)

$$\mathrm{Div}\left(\mathit{Q}\right)+J=0,$$

(5d)

where ${\mathit{Y}}_{\alpha }:=J{\hat{\mathit{y}}}_{\alpha }{\mathit{F}}^{-T}$ and $\mathit{Q}:=J\hat{\mathit{q}}{\mathit{F}}^{-T}$ denote, respectively, the backward Piola transform of the mass flux vector ${\mathit{y}}_{\alpha }$ and of the filtration velocity $\hat{\mathit{q}}:={\hat{\phi }}_{\mathrm{f}}[{\hat{\mathit{v}}}_{\mathrm{f}}-{\hat{\mathit{v}}}_{\mathrm{s}}]$ . We notice that equation (5b) can be obtained as a consequence of considering that the mass is conserved after growth has taken place [80]. Indeed, if we express both the mass of the tissue in its current configuration and in its natural state with respect to the reference configuration $\mathcal{B}$ , we will have that $J_{\mathrm{e}}{\hat{\phi }}_{\mathrm{s}}{\hat{\varrho }}_{\mathrm{s}}={\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_v$ , where ${\mathrm{\Phi }}_{\mathrm{s}v}:=J_{\mathrm{e}}{\hat{\phi }}_{\mathrm{s}}$ and ${\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_v$ denote, respectively, the volumetric fraction of the solid phase and the apparent mass density in the natural state. Moreover, from the the mass balance of the solid phase as a whole, see equation (4b), written in its material form, one may infer that $\stackrel{.}{\overline{J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_v}}=J{\hat{r}}_{\mathrm{s}}$ . Consequently, if the growth process preserves the density, the apparent mass density in the natural state is time-independent, which leads to equation (5b). We notice that these considerations can be reformulated by selecting ${\mathit{F}}_{\gamma }$ in a manner that ensures the relationship $\mathrm{tr}\left({\mathit{F}}_{\gamma }^{-1}{\stackrel{·}{\mathit{F}}}_{\gamma }\right)=J{\hat{r}}_{\mathrm{s}}/\left(J_{\gamma }{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}\right)$ to be satisfied, thereby representing a constraint on ${\mathit{F}}_{\gamma }$ [75, 86].

Prior to proceeding, it is worth mentioning that in different models of tumour growth, the assignment of $R_{\gamma }$ is consistently based on phenomenological considerations (see, for instance, [15, 17, 35, 87). Here, however, in line with DiCarlo and Quiligotti [36], Epstein and Maugin [79], Cermelli et al. [88], and Grillo et al. [89], we adopt the methodology put forward in Grillo and Di Stefano [70, 71]. The approach, which relies on the admission of the existence of non-conventional forces [36] to explain the chemo-mechanical coupling and other essential factors involved in growth, circumvents making any initial assumptions about $R_{\gamma }:=J{\hat{r}}_s/\left(J_{\gamma }{\rho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}\right)$ , called growth law in Grillo and Di Stefano [70, 71]. Instead, following Grillo and Di Stefano [70, 71], $R_{\gamma }$ is derived by solving the dynamics of the growing medium under consideration, eliminating the need for any preconceived notions. Within this framework, we remark that equations (5a)–(5d) can only be solved once the law for $R_{\gamma }$ has been found since $J{\hat{r}}_{\mathrm{s}}$ depends on $R_{\gamma }$ through the relationship $J{\hat{r}}_{\mathrm{s}}=J_{\gamma }{\rho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{R}_{\gamma }$ . Furthermore, by utilising the aforementioned relationship and the definition of r_s, once the expression for $R_{\gamma }$ is found, it will allow us to define the source/sink terms ${\hat{r}}_{\mathrm{fp}},{\hat{r}}_{\mathrm{nf}},{\hat{r}}_{\mathrm{ph}},\mathrm{and}{\hat{r}}_{\mathrm{nh}}$ . It is worth noting that unlike in Di Stefano et al. [24], the definitions of ${\hat{r}}_{\mathrm{ph}}\mathrm{and}{\hat{r}}_{\mathrm{nh}}$ are not provided a priori and will be established a posteriori based on the representation of non-conventional forces. Finally, we mention that a different viewpoint has been put forward in Grillo and Di Stefano [86] in which the growth law is treated as a non-holonomic kinematic constraint and reconsidered in Grillo and Di Stefano [90] without adopting time as a fictitious Lagrangian parameter.

3.2. Governing dynamic equations

Building upon the research conducted in DiCarlo and Quiligotti [36], Grillo and Di Stefano [70], and Grillo et al. [75], we adopt the perspective that the growth tensor serves as a kinematic variable that captures the tumour’s inelastic distortions induced by its growth. Although it is possible to account for all the components of ${\mathit{F}}_{\gamma }$ to deduce the equations that govern their dynamics as in DiCarlo and Quiligotti [36] and Grillo and Di Stefano [70]. In this work, we have assumed that the growth tensor ${\mathit{F}}_{\gamma }$ is spherical, namely, ${\mathit{F}}_{\gamma }=\gamma \mathit{I}$ , so that we seek to determine both the growth parameter γ and the growth law $R_{\gamma }$ by introducing suitable non-conventional forces [36, 70]. Within this context, we retrieve the equations determining the dynamics of the mixture by adhering to the Principle of Virtual Work and incorporate the set of kinematic descriptors {χ, F , γ} alongside their virtual variations, denoted by {δχ, δ F , δγ}. Then, for negligible inertial terms and no body forces, in the reference configuration, we can write [36, 70, 71]:

$$\begin{array}{c}\multicolumn{1}{c}{\sum _{\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}}\{{\int }_{\mathcal{B}}{\mathit{P}}_k:\mathit{g}\mathrm{Grad}{\delta }_{\chi }\mathrm{d}V+{\int }_{\mathcal{B}}J{\hat{\mathit{m}}}_{\mathrm{k}}.\delta \chi \mathrm{d}V\}+{\int }_{\mathcal{B}}\eta J_{\gamma }^{-1}\delta J_{\gamma }\mathrm{d}V}\\ \multicolumn{1}{c}{=\sum _{\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}}\{{\int }_{\partial \mathrm{N}\mathcal{B}}{\mathit{\tau }}_{\mathrm{k}}·{\delta }_{\chi }\mathrm{d}A\}+{\int }_{\mathcal{B}}\zeta J_{\gamma }^{J-1}\delta J_{\gamma }\mathrm{d}V,}\end{array}$$

(6)

where, for $\mathrm{k}\in \{\mathrm{s},\mathrm{f}\},{\mathit{P}}_{\mathrm{k}}:=J{\hat{\sigma }}_{\mathrm{k}}{\mathit{F}}^{-\mathrm{T}}$ is the first Piola–Kirchhoff stress tensor of the kth phase, g is the metric tensor associated with ${\mathcal{B}}_t$ , ${\hat{\mathit{m}}}_{\mathrm{k}}$ denote the internal forces describing the gain or loss of momentum of the kth phase in response to exchange interactions between the two phases, and τ _k is a boundary contact force on the Neumann portion of the boundary of $\mathcal{B}$ . Furthermore, within the framework of this work, the non-conventional forces η and ζ represent the spherical parts of the so-called remodelling couples [36, 70, 71] and are work conjugated of $J_{\gamma }^{-1}\delta J_{\gamma }$ . According to Di Carlo [91], it is worth noting that the external non-conventional force ζ has the potential to encompass interactions originating from smaller scales, such as biochemical interactions. Finally, we notice that our model is of order zero in the growth parameter γ, which describes the evolution of the inelastic distortions, while it is of order one in χ, which represents the observable deformation of the body.

In particular, equation (6) can be equivalently rewritten as:

$$\sum _{\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}}\left\{{\int }_{{\partial }_{\mathrm{N}}\mathcal{B}}[{\mathit{P}}_{\mathrm{k}}:\mathit{N}-{\mathit{\tau }}_{\mathrm{k}}]{\delta }_{\chi }\mathrm{d}A\right\}+\sum _{\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}}\{\underset{\mathcal{B}}{\int }[J{\hat{\mathit{m}}}_{\mathrm{k}}-\mathrm{Div}{\mathit{P}}_{\mathrm{k}}]·\delta \chi \mathrm{d}V\}+{\int }_{\mathcal{B}}[\eta -\zeta ]J_{\gamma }^{-1}\delta J_{\gamma }\mathrm{d}V=0,$$

(7)

so that, by localisation and adding a Dirichlet boundary condition for χ, we obtain the following system of equations:

$$J{\hat{\mathit{m}}}_{\mathrm{k}}-\mathrm{Div}{\mathit{P}}_{\mathrm{k}}=0,\mathrm{in}\mathcal{B},$$

(8a)

$${\mathit{P}}_{\mathrm{k}}·\mathit{N}={\mathit{\tau }}_{\mathrm{k}},\mathrm{on}{\partial }_{\mathrm{N}}\mathcal{B},$$

(8b)

$$\chi ={\chi }_{\mathrm{b}},\mathrm{on}{\partial }_{\mathrm{D}}\mathcal{B},$$

(8c)

$$\eta =\zeta ,\mathrm{in}\mathcal{B},$$

(8d)

for each $\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}$ and with ${\partial }_{\mathrm{D}}\mathcal{B}$ being the Dirichlet portion of the boundary of $\mathcal{B}$ . Therefore, by replacing equation (8a) for k = s with the equation obtained by summing (8a) over all $\mathrm{k}\in \{\mathrm{s},\mathrm{f}\}$ , we can express equations (8a)–(8d) as:

$$\mathrm{Div}\mathit{P}=0,\mathrm{in}\mathcal{B},$$

(9a)

$$J{\hat{\mathit{m}}}_{\mathrm{f}}-\mathrm{Div}{\mathit{P}}_{\mathrm{f}}=0,\mathrm{in}\mathcal{B},$$

(9b)

$$\mathit{P}.\mathit{N}=\tau ,\mathrm{on}{\partial }_{\mathrm{N}}\mathcal{B},$$

(9c)

$$\chi ={\chi }_{\mathrm{b}},\mathrm{on}{\partial }_{\mathrm{D}}\mathcal{B},$$

(9d)

$$\eta =\zeta ,\mathrm{in}\mathcal{B}.$$

(9e)

We notice that equation (9a) is the balance of linear momentum for the mixture as a whole, where $\mathit{P}:={\mathit{P}}_{\mathrm{s}}+{\mathit{P}}_{\mathrm{f}}$ is the first Piola-Kirchhoff stress tensor of the mixture and equation (9c) is the associated Neumann boundary condition with $\tau :={\tau }_{\mathrm{s}}+{\tau }_{\mathrm{f}}$ . In writing equation (9a), we made use of the condition $J({\hat{\mathit{m}}}_{\mathrm{s}}+{\hat{\mathit{m}}}_{\mathrm{f}})=0$ , which is a consequence of the principle of material frame-indifference [35, 76]. Furthermore, equation (9e), which represents a balance of growth-conjugated forces, has garnered substantial attention in numerous studies focused on growth and remodelling (see, for instance, [35, 92, 93]) and, by utilising the dissipation inequality as a method to establish a constitutive law governing the force η, it serves as the fundamental basis for deriving a dynamic equation that governs the behaviour of γ.

4. Study of the dissipation

In our analysis for obtaining constitutive information, we turn our attention to the compliance of P _s, P _f, η and Q with the dissipation inequality, which here we adapt from Grillo et al. [35] and Grillo and Di Stefano [70]. For clarity and convenience, we express the dissipation inequality in relation to the reference configuration as:

$$0\le {\int }_{\mathcal{C}}{\mathcal{D}}_{\mathrm{R}}\mathrm{d}V=-{\int }_{\mathit{C}}J{\hat{r}}_{\mathrm{s}}({\hat{\psi }}_{\mathrm{s}}-{\hat{\psi }}_{\mathrm{f}})\mathrm{d}V-{\int }_{\mathcal{C}}J{\varrho }_{\mathrm{s}0}{\hat{\phi }}_{\mathrm{s}}{\stackrel{·}{\hat{\psi }}}_{\mathrm{s}}+{\int }_{\mathcal{C}}{\mathrm{\Phi }}_{\mathrm{f}}^{-1}\mathit{Q}{\mathit{M}}_{\mathrm{f}}\mathrm{d}V+{\int }_{\mathcal{C}}{\mathit{P}}_{\mathrm{s}}:\mathit{g}\stackrel{·}{\mathit{F}}\mathrm{d}V+{\int }_{\mathcal{C}}{\mathit{P}}_{\mathrm{f}}:\mathit{g}\mathrm{Grad}{\hat{\mathit{v}}}_f\mathrm{d}V+{\int }_{\mathcal{C}}\zeta J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }\mathrm{d}V,$$

(10)

where $\mathcal{C}$ is an arbitrary open subset of $\mathcal{B}$ , ${\hat{\psi }}_{\mathrm{s}}$ and ${\hat{\psi }}_{\mathrm{f}}$ are the material form of the Helmholtz free energies per unit mass of the solid and fluid phases, and ${\mathit{M}}_{\mathrm{f}}:=J\left(\mathit{g}{\hat{\mathit{m}}}_{\mathrm{f}}\right)\mathit{F}$ is the material counterpart of the momentum exchange rate. In particular, as in Grillo et al. [75], in writing equation (10), we have approximated ${\hat{\psi }}_{\mathrm{f}}$ with a constant. Furthermore, the last term in the dissipation inequality represents the external net power conjugated with $J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }$ , which by means of equation (6), can be put in terms of the internal power [70, 71]. We further observe that equation (10) does not account for the dissipative contribution associated with the nutrients. Its omission is due to the complexities that would arise if one wants to derive a thermodynamically valid non-local constitutive law for the mass flux vector. Doing so is currently a subject of investigation in our research, but that falls outside the scope of this work.

By regarding the balance of mass of the mixture as a whole (i.e., the equation resulting from adding together equations (4b) and (4d)) as a constraint for the dissipation inequality with associated Lagrange multiplier p, the inequality (10) can be equivalently rewritten as:

$$\begin{array}{c}\multicolumn{1}{c}{0\le {\int }_{\mathcal{C}}{\mathcal{D}}_{\mathrm{R}}\mathrm{d}V=-{\int }_{\mathit{C}}J{\hat{r}}_{\mathrm{s}}({\hat{\psi }}_{\mathrm{s}}-{\hat{\psi }}_{\mathrm{f}})\mathrm{d}V-{\int }_{\mathcal{C}}J{\varrho }_{\mathrm{s}0}{\hat{\phi }}_{\mathrm{s}}{\stackrel{·}{\hat{\psi }}}_{\mathrm{s}}\mathrm{d}V+{\int }_{\mathcal{C}}{\mathrm{\Phi }}_{\mathrm{f}}^{-1}\mathit{Q}{\mathit{M}}_{\mathrm{f}}\mathrm{d}V+{\int }_{\mathcal{C}}{\mathit{P}}_{\mathrm{s}}:\mathit{g}\stackrel{·}{F}\mathrm{d}V+{\int }_{\mathcal{C}}{\mathit{P}}_{\mathrm{f}}:\mathit{g}\mathrm{Grad}{\hat{\mathit{v}}}_f\mathrm{d}V+{\int }_{\mathcal{C}}\zeta J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }\mathrm{d}V}\\ \multicolumn{1}{c}{+{\int }_{\mathcal{C}}\hat{p}[{\mathrm{\Phi }}_{\mathrm{s}}{\mathit{F}}^{-\mathrm{T}}:\stackrel{·}{\mathit{F}}+{\mathrm{\Phi }}_{\mathrm{f}}{\mathit{F}}^{-T}:\mathrm{Grad}{\hat{\mathit{v}}}_{\mathrm{f}}+J{\mathrm{\Phi }}_{\mathrm{f}}^{-1}\mathit{Q}\mathrm{Grad}\left(J^{-1}{\mathrm{\Phi }}_{\mathrm{f}}\right)]\mathrm{d}V,}\end{array}$$

(11)

where ${\mathrm{\Phi }}_{\mathrm{k}}:=J{\hat{\phi }}_{\mathrm{k}}$ . In particular, to account for the inelastic distortions induced by growth and under the hypothesis of isotropic and hyperelastic material, the strain energy density ${\hat{\psi }}_{\mathrm{s}}$ is constitutively expressed as a function of F and ${\mathit{F}}_{\gamma }$ such that ${\hat{\psi }}_{\mathrm{s}}={\psi }_v^{⌣}\left({\mathit{FF}}_{\gamma }^{-1}\right)={\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)$ (see, for instance, [70, 75]). Thus, by means of the chain rule and taking into account (4b), the localisation of the dissipation inequality leads to:

$$\begin{array}{c}\multicolumn{1}{c}{0\le {\mathcal{D}}_{\mathrm{R}}={\mathrm{\Phi }}_{\mathrm{f}}^{-1}\mathit{Q}[{\mathit{M}}_{\mathrm{f}}+J\hat{p}\mathrm{Grad}\left(J^{-1}{\mathrm{\Phi }}_{\mathrm{f}}\right)]+[{\mathit{P}}_{\mathrm{f}}+\hat{p}{\mathrm{\Phi }}_{\mathrm{f}}{\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}}]:\mathit{g}\mathrm{Grad}{\hat{\mathit{v}}}_{\mathrm{f}}}\\ \multicolumn{1}{c}{+[{\mathit{P}}_{\mathrm{s}}+\hat{p}{\mathrm{\Phi }}_{\mathrm{s}}{\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}}-J_{\gamma }^{2/3}{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{\mathit{g}}^{-1}\frac{\partial {\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)}{\partial {\mathit{F}}_{\mathrm{e}}}{\mathit{I}}^{-\mathrm{T}}]:\mathit{g}\stackrel{·}{F}}\\ \multicolumn{1}{c}{+[\eta -\frac{1}{3}{J}_{\gamma }{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}(({\hat{\psi }}_{\mathrm{s}}-{\hat{\psi }}_{\mathrm{f}}){\mathit{I}}^{-\mathrm{T}}-J_{\gamma }^{-1/3}{\mathit{F}}^{\mathrm{T}}\frac{\partial {\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)}{\partial {\mathit{F}}_{\mathrm{e}}}{\mathit{I}}^{-\mathrm{T}}):\mathit{I}]J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma },}\end{array}$$

(12)

which by means of the Coleman and Noll procedure lead to the following identifications for the first Piola– Kirchhoff stress tensor of the solid and of the fluid phase:

$${\mathit{P}}_{\mathrm{s}}=-\hat{p}{\mathrm{\Phi }}_{\mathrm{s}}{\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}}+{\mathit{P}}_{\mathrm{sc}},$$

(13a)

$${\mathit{P}}_{\mathrm{f}}=-\hat{p}{\mathrm{\Phi }}_{\mathrm{f}}{\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}},$$

(13b)

where ${\mathit{P}}_{\mathrm{sc}}:=J_{\gamma }^{2/3}{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{\mathit{g}}^{-1}\frac{\partial {\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)}{\partial {\mathit{F}}_{\mathrm{e}}}{\mathit{I}}^{-\mathrm{T}}$ represents the constitutive part of the first Piola–Kirchhoff stress tensor of the solid phase [75]. Furthermore, by introducing the notations ${\mathit{M}}_{\mathrm{fd}}$ and ${\eta }_{\mathrm{d}}$ for the dissipative parts of the momentum exchange rate and the internal conjugated force, namely:

$${\mathit{M}}_{\mathrm{fd}}:={\mathit{M}}_{\mathrm{f}}+J\hat{p}\mathrm{Grad}\left(J^{-1}{\mathrm{\Phi }}_{\mathrm{f}}\right),$$

(14a)

$${\eta }_{\mathrm{d}}:=\eta -\mathcal{H},$$

(14b)

with:

$$\mathcal{H}:=J_{\gamma }{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}({\hat{\psi }}_{\mathrm{s}}-{\hat{\psi }}_{\mathrm{f}})-\frac{1}{3}\mathrm{tr}\left({\mathit{F}}^{\mathrm{T}}\mathit{g}{\mathit{P}}_{\mathrm{sc}}\right),$$

(15)

the local dissipation inequality rewrites as:

$$0\le {\mathcal{D}}_{\mathrm{R}}={\eta }_{\mathrm{d}}{J}_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }+{\mathrm{\Phi }}_{\mathrm{f}}^{-1}\mathit{Q}{\mathit{M}}_{\mathrm{fd}}.$$

(16)

We notice that equation (15) represents, up to a constant, the spherical part of the Eshelby stress tensor [70, 79, 94, 95, a tensor that has been identified to be the driving force of inelastic processes contributing to growth [79].

5. Constitutive considerations

Here, we set the constitutive considerations of our model, which involve formulating relationships between different quantities, such as fluxes and forces, that govern the system’s dynamics.

5.1. Darcy’s law

By adhering to established practices in the field of porous media mechanics (see, for example, [30, 35, 96]), we express ${\mathit{M}}_{\mathrm{fd}}$ constitutively as a linear function of the fluid filtration velocity, namely:

$${\mathit{M}}_{\mathrm{fd}}:={\mathrm{\Phi }}_{\mathrm{f}}{\mathit{K}}^{-1}\mathit{Q},$$

(17)

where $\mathit{K}:=J{\mathit{F}}^{-1}\hat{\mathit{k}}{\mathit{F}}^{-\mathrm{T}}$ , with $\hat{k}$ being the material form of the spatial permeability tensor, denotes the material permeability tensor, which is assumed to be invertible, symmetric, and positive definite. Therefore, the substitution of equations (13b) and (17) in equation (9b) yields:

$${\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}}[{\mathrm{\Phi }}_{\mathrm{f}}{\mathit{K}}^{-1}\mathit{Q}-J\hat{p}\mathrm{Grad}\left(J^{-1}{\mathrm{\Phi }}_{\mathrm{f}}\right)]-\mathrm{Div}{\mathit{P}}_{\mathrm{f}}=0,$$

(18)

and, consequently, using the Piola identity, equation (18) reduces to:

$$\mathit{Q}=-\mathit{K}\mathrm{Grad}\hat{p},$$

(19)

which represents the material form of Darcy’s law of filtration. As in previous works [24, 26], the permeability tensor of the tissue is assumed to be proportional to the inverse of the metric tensor, so that $\hat{\mathit{k}}=k_0{\mathit{g}}^{-1}$ . Then, we can write:

$$\mathit{K}=J{\mathit{F}}^{-1}{k}_0{\mathit{g}}^{-1}{\mathit{F}}^{-\mathrm{T}}=J{k}_0{\mathit{C}}^{-1},$$

(20)

where C is the right Cauchy-Green deformation tensor and:

$$k_0=k_R{\left(\frac{J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{sv}}}{J_{\gamma }{\phi }_{\mathrm{fR}}}\right)}^{m_0}\exp \left(\frac{m_1}{2}\left[\frac{J^2-J_{\gamma }^2}{J_{\gamma }^2}\right]\right),$$

(21)

is the Holmes–Mow scalar permeability [97]. In equation (21), $k_{\mathrm{R}}$ is a reference permeability, ${\phi }_{\mathrm{fR}}:=1-{\mathrm{\Phi }}_{\mathrm{s}v}$ denotes a reference value of the fluid phase volumetric fraction, and $m_0$ and $m_1$ are constant material parameters. Finally, the substitution of equations (20) and (21) into equation (19) leads to:

$$\mathit{Q}=-J\left\{k_{\mathrm{R}}{\left(\frac{J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}}{J_{\gamma }{\phi }_{\mathrm{fR}}}\right)}^{m_0}\exp \left(\frac{m_1}{2}\left[\frac{J^2-J_{\gamma }^2}{J_{\gamma }^2}\right]\right)\right\}{\mathit{C}}^{-1}\mathrm{Grad}\hat{p}.$$

(22)

5.2. Constitutive part of the first Piola–Kirchhoff stress tensor

Following Di Stefano et al. [24] and Ramírez-Torres et al. [26], the solid phase Helmholtz free energy density is assumed to be of the Holmes–Mow type [97], namely,

$${\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)={\overline{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right)={\tilde{\psi }}_v\left({\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}\right)=\frac{a_0}{{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}}[\exp \left(\mathrm{\Psi }\left({\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}\right)\right)-1],$$

(23a)

$$\mathrm{\Psi }\left({\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}\right)=a_1[I_1\left({\mathit{C}}_{\mathrm{e}}\right)-3]+a_2[I_2\left({\mathit{C}}_{\mathrm{e}}\right)-3]-a_3\log \left(I_3\left({\mathit{C}}_{\mathrm{e}}\right)\right),$$

(23b)

where ${\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}$ denotes the list of the principal invariants of the elastic right Cauchy–Green deformation tensor ${\mathit{C}}_{\mathrm{e}}:={\mathit{F}}_{\mathrm{e}}^{\mathrm{T}}\mathit{g}{\mathit{F}}_{\mathrm{e}}={\mathit{F}}_{\gamma }^{-\mathrm{T}}{\mathit{CF}}_{\gamma }^{-\mathrm{T}}$ , namely, ${\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}=(I_1\left({\mathit{C}}_e\right),I_2\left({\mathit{C}}_e\right),I_3\left({\mathit{C}}_e\right))$ , and a₀, a₁, a₂, and a₃ are defined as:

$$a_0:=\frac{2\mu +\lambda }{4a_3},a_1:=a_3\frac{2\mu -\lambda }{2\mu +\lambda },a_2:=a_3\frac{\lambda }{2\mu +\lambda },a_3:=a_1+2a_2=1$$

(24)

with λ and μ being Lamé’s parameters. The above considerations imply that the constitutive part of the first Piola–Kirchhoff stress tensor is given by:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{P}}_{\mathrm{sc}}=J_{\gamma }^{2/3}{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{\mathit{g}}^{-1}\frac{\partial {\psi }_v\left({\mathit{F}}_{\mathrm{e}}\right)}{\partial {\mathit{F}}_{\mathrm{e}}}{\mathit{I}}^{-\mathrm{T}}=2J_{\gamma }^{2/3}{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{\mathit{F}}_{\mathrm{e}}\frac{\partial {\psi }_v\left({\mathit{C}}_{\mathrm{e}}\right)}{\partial {\mathit{C}}_{\mathrm{e}}}{\mathit{I}}^{-\mathrm{T}}}\\ \multicolumn{1}{c}{=2J_{\gamma }\{b_1{\mathit{FC}}_{\gamma }^{-1}+b_2[I_1\left({\mathit{C}}_{\mathrm{e}}\right){\mathit{FC}}_{\gamma }^{-1}-{\mathit{FC}}_{\gamma }^{-1}{\mathit{CC}}_{\gamma }^{-1}]+b_3{I}_3\left({\mathit{C}}_{\mathrm{e}}\right)\mathit{F}{\mathit{C}}^{-1}\},}\end{array}$$

(25)

where:

$$b_1:=\frac{\partial {\tilde{\psi }}_v\left({\mathcal{I}}_{{\mathit{C}}_e}\right)}{\partial I_1\left({\mathit{C}}_{\mathrm{e}}\right)}=a_1[{\stackrel{⌣}{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right)+a_0],$$

(26a)

$$b_2:=\frac{\partial {\tilde{\psi }}_v\left({\mathcal{I}}_{{\mathit{C}}_e}\right)}{\partial I_2\left({\mathit{C}}_{\mathrm{e}}\right)}=a_2[{\stackrel{⌣}{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right)+a_0],$$

(26b)

$$b_3:=\frac{\partial {\tilde{\psi }}_v\left({\mathcal{I}}_{{\mathit{C}}_e}\right)}{\partial I_3\left({\mathit{C}}_{\mathrm{e}}\right)}=a_3{\left[I_3\left({\mathit{C}}_{\mathrm{e}}\right)\right]}^{-1}[{\stackrel{⌣}{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right)+a_0],$$

(26c)

where ${\stackrel{⌣}{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right):={\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{\overline{\psi }}_v\left({\mathit{C}}_{\mathrm{e}}\right)=a_0[\exp \left(\mathrm{\Psi }\left({\mathcal{I}}_{{\mathit{C}}_{\mathrm{e}}}\right)\right)-1].$

5.3. Growth law

Drawing inspiration from Grillo and Di Stefano [70], here, we recast an equation for the evolution of the growth parameter γ from the balance of growth-conjugated forces specified in equation (9e). To accomplish this, we begin by introducing the dissipative part of η in the natural state as ${\eta }_{\mathrm{d}v}=J_{\gamma }^{-1}{\eta }_{\mathrm{d}}$ , which, by referring to equations (14b) and (9e), leads to the expression:

$$\eta \mathrm{d}v=J_{\gamma }^{-1}[-\mathcal{H}+\zeta ],$$

(27)

where $\mathcal{H}$ has been defined in equation (15). Thus, to comply with the dissipation inequality, we set:

$${\eta }_{\mathrm{d}v}:=\frac{1}{3}{a}_v{J}_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma },$$

(28)

where a_v is a positive real constant with physical dimensions of stress per unit of time. Consequently, the sought dynamic equation for γ can be written in the form:

$$\frac{1}{3}{a}_v{J}_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }=J_{\gamma }^{-1}[-\mathcal{H}+\zeta ],$$

(29)

or equivalently, using equation (3), as:

$${\gamma }^{-1}\stackrel{·}{\gamma }=a_v^{-1}{J}_{\gamma }^{-1}[-\mathcal{H}+\zeta ].$$

(30)

We remark that, as shown in Grillo and Di Stefano [70], a law for ${\eta }_{\mathrm{d}v}$ that satisfies fundamental principles of constitutive laws [98], such as independence of the choice of the reference configuration, can be written in the form given in equation (28). Furthermore, by examining equation (30), we observe that the dynamics of γ are determined by the difference between $\mathcal{H}$ and ζ. Consequently, even when $\zeta \equiv 0$ , the governing equation for γ remains non-trivial due to the influence of $\mathcal{H}$ , which drives the evolution of the growth parameter. Nevertheless, when a non-zero ζ is present, its specific form plays a crucial role in determining how γ evolves. Among the possible choices for ζ and, since we are dealing with the growth of an avascular tumour where different constituents are interacting, we assign an expression accounting for the nutrients and mechanotransduction. In particular, following a methodology similar to Grillo and Di Stefano [70] and by targeting the phenomenological expressions found in Mascheroni et al. [9], Di Stefano et al. [24], and Ramírez-Torres et al. [26], so that we preserve the biological concepts and ideas supported by previous research efforts and that have been successful in the modelling of tumour growth, we write:

$$\zeta :=J_{\gamma }({\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}),$$

(31)

where ${\kappa }_{\mathrm{fp}}$ and ${\kappa }_{\mathrm{nf}}$ are the constant parameters with physical units of stress per time, and the scalar quantities ${\mathrm{\Upsilon }}_{\mathrm{fp}}$ and ${\mathrm{\Upsilon }}_{\mathrm{nf}}$ are defined as:

$${\mathrm{\Upsilon }}_{\mathrm{fp}}:=\frac{{\zeta }_{\mathrm{fp}}}{3{\varrho }_{\mathrm{s}0}}{〈\frac{{\hat{\omega }}_{\mathrm{cf}}-{\omega }_{\mathrm{cr}}}{{\omega }_{\mathrm{env}}-{\omega }_{\mathrm{cr}}}〉}_{+}[1-\frac{{\delta }_1{〈\overline{\sigma }〉}_{+}}{{\delta }_2+{〈\overline{\sigma }〉}_{+}}]\frac{J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}}{J{\phi }_{\mathrm{fR}}}{\hat{\omega }}_{\mathrm{ps}},$$

(32a)

$${\mathrm{\Upsilon }}_{\mathrm{nf}}:=\frac{{\zeta }_{\mathrm{nf}}}{3{\varrho }_{\mathrm{s}0}}(1-{\hat{\omega }}_{\mathrm{ps}}).$$

(32b)

In equations (32a) and (32b), ${\zeta }_{\mathrm{fp}}$ and ${\zeta }_{\mathrm{nf}}$ account for the absorption of fluid by the proliferating cells and for the necrotic cells that dissolve into the fluid, ${\omega }_{\mathrm{cr}}$ is a critical value for the nutrients’ mass fraction below which proliferating cells start dying and ${\omega }_{\mathrm{env}}$ represents the value of the nutrients available in the tissue’s environment. Furthermore, ${\delta }_1$ is a dimensionless constant, ${\delta }_2$ is a strictly positive characteristic stress, and $\overline{\sigma }$ denotes the spherical part of the Cauchy stress tensor. Finally, the notation ${〈\varphi 〉}_{+}$ is used to represent Macaulay’s bracket that returns zero if $\varphi$ is negative or zero, and $\varphi$ if it is positive.

Thus, combining the above results, the evolution equation for the growth parameter, γ, is:

$${\gamma }^{-1}\stackrel{·}{\gamma }=-\frac{\mathcal{H}}{a_v{J}_{\gamma }}+\frac{1}{a_v}({\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}).$$

(33)

Furthermore, from equations (5b) and (30), we are now able to provide an expression for the growth law $R_{\gamma }$ , which is similar to the one found in Grillo and Di Stefano [70] and reads:

$$R_{\gamma }\equiv J_{\gamma }^{-1}{\stackrel{·}{J}}_{\gamma }=-\frac{3\mathcal{H}}{a_v{J}_{\gamma }}+\frac{3}{a_v}({\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}).$$

(34)

So, once the growth parameter and the other unknowns of the problem have been identified, the value of $R_{\gamma }$ can be determined through equation (34). In addition, by recalling that $J{\hat{r}}_{\mathrm{s}}=J_{\gamma }{\varrho }_{\mathrm{s}0}{\mathrm{\Phi }}_{\mathrm{s}v}{R}_{\gamma }$ with $J{\hat{r}}_{\mathrm{s}}=J[{\hat{r}}_{\mathrm{fp}}+{\hat{r}}_{\mathrm{nf}}+{\hat{r}}_{\mathrm{ph}}+{\hat{r}}_{\mathrm{nh}}]$ , we can make the following identifications:

$$\frac{J{\hat{r}}_{\mathrm{fp}}}{J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}}\equiv \frac{3{\kappa }_{\mathrm{fp}}}{a_v}{\mathrm{\Upsilon }}_{\mathrm{fp}}=\frac{{\zeta }_{\mathrm{fp}}{\kappa }_{\mathrm{fp}}}{{\varrho }_{\mathrm{s}0}{a}_v}{〈\frac{{\hat{\omega }}_{\mathrm{cf}}-{\omega }_{\mathrm{cr}}}{{\omega }_{\mathrm{env}}-{\omega }_{\mathrm{cr}}}〉}_{+}[1-\frac{{\delta }_1〈\overline{\sigma }〉}{{\delta }_2+{〈\overline{\sigma }〉}_{+}}]\frac{J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}}{J{\phi }_{\mathrm{fr}}}{\hat{\omega }}_{\mathrm{ps}},$$

(35a)

$$\frac{J{\hat{r}}_{\mathrm{nf}}}{J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}}\equiv \frac{3{\kappa }_{\mathrm{nf}}}{a_v}{\mathrm{\Upsilon }}_{\mathrm{nf}}=-\frac{{\zeta }_{\mathrm{nf}}{\kappa }_{\mathrm{nf}}}{{\varrho }_{\mathrm{s}0}{a}_v}(1-{\hat{\omega }}_{\mathrm{ps}}),$$

(35b)

$$\frac{J{\hat{r}}_{\mathrm{ph}}}{J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}}=\frac{3\mathcal{H}}{a_v{J}_{\gamma }},$$

(35c)

where we have set ${\hat{r}}_{\mathrm{nh}}\equiv 0$ . As remarked in previous sections, equation (35c) establishes an a posteriori connection between the production of inelastic distortions, encoded in the growth parameter γ, and ${\hat{r}}_{\mathrm{ph}}$ , which reflects, by referring to equation (5a), the impact of the properties of ${\mathit{F}}_{\gamma }$ on the evolution of the proliferating cells.

5.4. Mass flux vector

By taking inspiration from Ramírez-Torres et al. [26], the components of the mass flux vector of nutrients are defined through the expression:

$${\left[{\mathit{y}}_{\alpha }(x,t)\right]}^a=-\varrho f{\int }_{{\mathcal{B}}_t}{f}_{\alpha (x,\tilde{x})}{\left[\mathit{d}(x,\tilde{x},t)\right]}^{\mathrm{ab}}{\left[\mathrm{grad}\left({\omega }_{\mathrm{cf}}(\tilde{x},t)\right)\right]}_b\mathrm{d}v\left(\tilde{x}\right),$$

(36)

where $f_{\alpha (x,\tilde{x})}(x-\tilde{x})$ is a non-locality function depending on the dimensionless positive scalar function $\alpha (x,\tilde{x})\in {\mathbb{R}}^{+}$ , and $\mathit{d}(x,\tilde{x},t)$ is the diffusivity tensor. Moreover, we notice that since ${\mathit{y}}_{\alpha }$ represents the flux of nutrients, it needs to have physical units of a standard flux, i.e., $M{L}^{-2}{T}^{-1}$ , where M, L, and T denote a specific unit of mass, length, and time, respectively. This implies that the non-locality function has physical units L^–3. We remark that the non-uniform form of the scalar quantity $\alpha$ represents an additional feature with respect to the definition given in Ramírez-Torres et al. [26].

Remark 1 By considering the notions offered by fractional calculus [48–50], we observe that it is possible to assume the non-locality function to be a power-law dependent on the distance between the two spatial points x and $\tilde{x}$ , with $(x,\tilde{x})$ being the exponent, such that the integral in equation (36) can be thought as a variable-order fractional differentiation operator. In particular, within the context of anomalous diffusion, the consideration of a space-dependent order $\alpha (x,\tilde{x})$ has been linked to the medium heterogeneity (see, for instance, [61] and references therein) such that it can be interpreted as a quantity encoding location-dependent, microstructural properties of a system at larger scales.

In addition, we notice that, just like constant-order fractional operators, the variable-order counterparts lack a universally accepted definition, resulting in numerous options presented in the scientific literature [59, 64, 99–105]. The suggested alternatives reflect desired mathematical properties or physical interpretations under specific conditions. Emphasis has been placed on the functional expression of a time-dependent fractional differentiation order, aiming to capture diverse memory patterns in the modelling process. Notably, with abuse of notation, the cases $\alpha (t,\tau )=\alpha \left(t\right)$ and $\alpha (t,\tau )=\alpha \left(t\right)$ offer different forms on how the fractional operator “remembers” past values of $\alpha$ , with the memory being weaker when $\alpha (t,\tau )=\alpha \left(t\right)$ , as it only depends on the current time, and stronger when $\alpha (t,\tau )=\alpha \left(t\right)$ as it contains the memory of the order history [60, 61]. Within the context of this work, however, we put emphasis on the spatial dependence in the variable-order to account for potential spatial phenomena that may arise due to the heterogeneous and complex environment in which diffusion is taking place. Thus, the aforementioned interpretations translate to accounting for locationdependent, non-local interactions determined by the current position when $\alpha (x,\tilde{x})=\alpha \left(x\right)$ or relying on the information at “distant” points when $\alpha (x,\tilde{x})=\alpha \left(x\right)$ . It is important to note that our system evolves over time, so the variable-order parameter $\alpha (x,\tilde{x})$ is influenced by the motion since $x=\chi (X,t)$ and $\tilde{x}=\chi (\tilde{X},t)$ . This reflects a type of variable-order non-locality that changes with the dynamics of the medium being studied.

From the definition of the material mass flux vector, we have that:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{Y}}_{\alpha }(X,t):=J(X,t){\hat{\mathit{y}}}_{\alpha }(X,t){\mathit{F}}^{-T}(X,t)}\\ \multicolumn{1}{c}{=-{\varrho }_f{\int }_{\mathcal{B}}{\Im }_{\alpha }(\chi (X,t)-\chi (\tilde{X},t))\mathfrak{D}(X,\tilde{X},t)\mathrm{Grad}\left({\hat{\omega }}_{\mathrm{cf}}(\tilde{X},t)\right)\mathrm{d}V\left(\tilde{X}\right),}\end{array}$$

(37)

where ${\Im }_{\alpha }$ is the material non-locality function and $\mathfrak{D}$ denotes the non-local material diffusivity tensor and are defined as:

$${\Im }_{\alpha (\chi (X,t),\chi (\tilde{X},t))}(\chi (X,t)-\chi (\tilde{X},t)):=J(X,t)f_{\alpha (\chi (X,t),\chi (\tilde{X},t))}(\chi (X,t)-\chi (\tilde{X},t)),$$

(38a)

$$\mathfrak{D}(X,\tilde{X},t):=J(\tilde{X},t){\mathit{F}}^{-1}(\chi (X,t),t)\hat{\mathit{d}}(X,\tilde{X},t){\mathit{F}}^{-\mathrm{T}}(\tilde{X},t).$$

(38b)

We observe that the material non-locality function considers a form of non-locality that changes according to the tumour’s kinematics through its motion χ.

Building on previous research [26], we suppose that the diffusivity tensor d depends on $\tilde{x}$ and t, so that if $d(\tilde{x},t)={\phi }_f(\tilde{x},t)d_{\mathrm{R}}{\mathit{g}}^{-1}(\tilde{x},t)$ , where ${\mathit{d}}_{\mathrm{R}}$ is a constant reference diffusivity [24, 26], we can write:

$$\hat{\mathit{d}}(\tilde{X},t)=\frac{J(\tilde{X},t)-J_{\gamma }(\tilde{X},t){\mathrm{\Phi }}_{\mathrm{sv}}}{J(\tilde{X},t)}{d}_{\mathrm{R}}{\mathit{g}}^{-1}(\chi (\tilde{X},t),t),$$

(39)

and the material diffusivity tensor given in equation (38b) can be expressed in the form:

$$\mathfrak{D}(X,\tilde{X},t)=[J(\tilde{X},t)-J_{\gamma }(\tilde{X},t){\mathrm{\Phi }}_{\mathrm{s}v}]d_{\mathrm{R}}{\mathit{F}}^{-1}(\chi (X,t),t){\mathit{g}}^{-1}(\chi (\tilde{X},t),t){\mathit{F}}^{-\mathrm{T}}(\tilde{X},t),$$

(40)

which reduces to the standard material diffusivity tensor, i.e., $\mathfrak{D}=[J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}]d_{\mathrm{R}}\mathit{C}$ , in the case in which X and $\tilde{X}$ are indistinguishable [26].

6. Summary of the model and benchmark problem

In this section, we provide an overview of the equations under consideration. In summary, we have:

$${\stackrel{·}{\hat{\omega }}}_{\mathrm{ps}}=J[{\hat{r}}_{\mathrm{pn}}+{\hat{r}}_{\mathrm{fp}}+{\hat{r}}_{\mathrm{ph}}]{\left[J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}\right]}^{-1}-R_{\gamma }{\hat{\omega }}_{\mathrm{ps}},$$

(41a)

$$\stackrel{·}{\gamma }=-\frac{\mathcal{H}}{a_v{\gamma }^2}+\frac{\gamma }{a_v}[{\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}],$$

(41b)

$$R_{\gamma }=-\frac{3\mathcal{H}}{a_v{J}_{\gamma }}+\frac{3}{a_v}[{\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}],$$

(41c)

$$[J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}]{\stackrel{·}{\hat{\omega }}}_{\mathrm{cf}}-\left[\mathit{K}\mathrm{Grad}\hat{p}\right]\mathrm{Grad}\left({\hat{\omega }}_{\mathrm{cf}}\right)+\frac{1}{{\varrho }_{\mathrm{s}0}}\mathrm{Div}\left({\mathit{Y}}_{\alpha }\right)=\frac{J{\hat{r}}_{\mathrm{cf}}}{{\varrho }_{\mathrm{s}0}}+J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{R}_{\gamma }{\hat{\omega }}_{\mathrm{cf}},$$

(41d)

$$\mathrm{Div}(-J\hat{p}{\mathit{g}}^{-1}{\mathit{F}}^{-T}+{\mathit{P}}_{\mathrm{sc}})=0,$$

(41e)

$$\stackrel{·}{J}-\mathrm{Div}\left(\mathit{K}\mathrm{Grad}\hat{p}\right)=0,$$

(41f)

where K , ${\mathit{P}}_{\mathrm{sc}}$ , $\mathcal{H}$ , and ${\mathit{Y}}_{\alpha }$ have been specified in equations (20), (25), (15), and (37), respectively, and as in Di Stefano et al. [24] and Ramírez-Torres et al. [26], we set:

$$\frac{J{\hat{r}}_{\mathrm{pn}}}{J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}}=-\frac{{\zeta }_{\mathrm{pn}}}{{\varrho }_{\mathrm{s}0}}{〈1-\frac{{\hat{\omega }}_{\mathrm{cf}}}{{\omega }_{\mathrm{cr}}}〉}_{+}{\hat{\omega }}_{\mathrm{ps}},$$

(42a)

$$\frac{J{\hat{r}}_{\mathrm{cf}}}{{\varrho }_{\mathrm{s}0}}=-\frac{{\zeta }_{\mathrm{cf}}}{{\varrho }_{\mathrm{s}0}}\frac{{\hat{\omega }}_{\mathrm{cf}}}{{\hat{\omega }}_{\mathrm{cf}}+{\omega }_0}{J}_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\hat{\omega }}_{\mathrm{ps}},$$

(42b)

with ${\omega }_0\in ]0,1]$ being a reference value influencing the consumption of nutrients, and ${\zeta }_{\mathrm{pn}}$ and ${\zeta }_{\mathrm{cf}}$ representing factors that contribute to the consumption of nutrients and the emergence of necrotic cells.

One significant distinction from the model presented in Ramírez-Torres et al. [26], is that for growth to occur, it is not necessary that $J{\hat{r}}_{\mathrm{fp}}$ is positive as it was the case in Ramírez-Torres et al. [26]. In the present framework, growth can occur even when $J{\hat{r}}_{\mathrm{fp}}\le 0$ as the contribution of the configurational mechanical stress $\mathcal{H}$ would still be present and driving the evolution of the body [70, 86].

6.1. Description of a benchmark problem

Now, we focus our attention on the description of the benchmark problem previously investigated in [24, 26, 106]. Specifically, we examine a scenario in which the growth of an avascular tumour occurs within a cylindrical domain with a rigid curved boundary, so that the tumour’s motion is exclusively constrained to happen in the axial direction. This description is aligned with the progression of DCIS, which represents the earliest stage of breast cancer. Indeed, DCIS is characterised by a non-invasive state in which tumour cells are confined by an intact basement membrane acting as a thin, flexible sheet-like extracellular matrix structure, effectively preventing the penetration of blood vessels from the surrounding tissue into the duct [73].

As in Di Stefano et al. [24] and Ramírez-Torres et al. [26], we describe the reference configuration of the tissue using a cylindrical coordinate system (R, Θ, Z), where $R\in [0,R_{\mathrm{in}}[$ , $\mathrm{\Theta }\in [0,2\pi [,$ and $Z\in [-L_{\mathrm{in}}+L_{\mathrm{in}}]$ represent the radial, circumferential, and axial coordinates, respectively. In particular, R_in denotes the initial radius and 2L_in is the initial length of the cylinder. Correspondingly, the current configuration of the tissue is characterized by the cylindrical coordinate system $(r,\theta ,z)$ , where:

$$r={\chi }^r(R,\mathrm{\Theta },Z,t)=R,$$

(43a)

$$\theta ={\chi }^{\theta }(R,\mathrm{\Theta },Z,t)=\mathrm{\Theta },$$

(43b)

$$z={\chi }^z(R,\mathrm{\Theta },Z,t)=Z+u(Z,t),$$

(43c)

with $u(Z,t)$ being the axial displacement. Specifically, we focus on maintaining a constant radius for the specimen, but allowing its length to vary solely along the axial direction. This approach effectively eliminates any rotational movement around the primary axis. Accordingly, the set of unknowns of the problem at hand, namely, $\mathcal{U}=\{{\hat{\omega }}_{\mathrm{ps}},\gamma ,{\hat{\omega }}_{\mathrm{cf}},u,\hat{p}\}$ depends only on the axial coordinate and time. Furthermore, the restrictions imposed to the motion imply that while equations (41a), (41b), and (41c) will have the same structure, equations (41d)–(41f) can be further simplified. In particular, in the framework of the benchmark problem, the system of equations given in the previous section reduces to:

$${\stackrel{·}{\hat{\omega }}}_{\mathrm{ps}}=J[{\hat{r}}_{\mathrm{pn}}+{\hat{r}}_{\mathrm{fp}}+{\hat{r}}_{\mathrm{ph}}]{\left[J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{\varrho }_{\mathrm{s}0}\right]}^{-1}-R_{\gamma }{\hat{\omega }}_{\mathrm{ps}},$$

(44a)

$$\stackrel{·}{\gamma }=-\frac{\mathcal{H}}{a_v{\gamma }^2}+\frac{\gamma }{a_v}[{\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}],$$

(44b)

$$R_{\gamma }=-\frac{3\mathcal{H}}{a_v{J}_{\gamma }}+\frac{3}{a_v}[{\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}],$$

(44c)

$$[J-J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}]{\stackrel{·}{\hat{\omega }}}_{\mathrm{cf}}-k_0{[1+\frac{\partial u}{\partial Z}]}^{-1}\frac{\partial \hat{p}}{\partial Z}\frac{\partial {\hat{\omega }}_{\mathrm{cf}}}{\partial Z}+\frac{1}{{\varrho }_{\mathrm{s}0}}\frac{\partial Y_{\alpha }^Z}{\partial Z}=\frac{J{\hat{r}}_{\mathrm{cf}}}{{\varrho }_{\mathrm{s}0}}+J_{\gamma }{\mathrm{\Phi }}_{\mathrm{s}v}{R}_{\gamma }{\hat{\omega }}_{\mathrm{cf}},$$

(44d)

$$\frac{\partial }{\partial Z}(-\hat{p}+(1+\frac{\partial u}{\partial Z}){\left[{\mathit{S}}_{\mathrm{sc}}\right]}^{\mathrm{ZZ}})=0,$$

(44e)

$$\frac{{\partial }^{2}u}{\partial t\partial Z}-\frac{\partial }{\partial Z}\left(k_0{(1+\frac{\partial u}{\partial Z})}^{-1}\frac{\partial \hat{p}}{\partial Z}\right)=0,$$

(44f)

where we have used the fact that $\left[\mathit{F}\right]=\mathrm{diag}\{1,1,1+\partial u/\partial Z\}$ and ${\mathit{S}}_{\mathrm{sc}}={\mathit{F}}^{-1}{\mathit{P}}_{\mathrm{sc}}$ is the second Piola–Kirchhoff stress tensor.

To solve the system of equations specified by equations (44a)–(44f), we use the same initial and boundary conditions as those employed in Di Stefano et al. [24] and Ramírez-Torres et al. [26]. Specifically, at the initial time instant $t=0$ , we assume that the solid phase consists entirely of proliferating cells, i.e., ${\hat{\omega }}_{\mathrm{ps}}(Z,0)=1$ , and that the growth parameter is uniform and equal to 1. Besides, nutrients are considered to be evenly distributed in the spatial domain with mass fraction ${\hat{\omega }}_{\mathrm{cf}}(Z,0)={\omega }_{\mathrm{env}}>0$ . Finally, both the pressure and the axial displacement are assumed to be zero at the initial time.

On the contrary, by denoting with ${(\partial \mathcal{B})}_{\mathrm{left}}$ and ${(\partial \mathcal{B})}_{\mathrm{right}}$ the left and right surfaces at the extremities of $\mathcal{B}$ , respectively, and with ${(\partial \mathcal{B})}_{\mathrm{C}}$ the lateral boundary of the cylinder, the boundary conditions complementing equations (44a)–(44f) are:

$${\hat{\omega }}_{\mathrm{cf}}={\omega }_{\mathrm{env}},\mathrm{on}{(\partial \mathcal{B})}_{\mathrm{left}}\mathrm{and}{(\partial \mathcal{B})}_{\mathrm{right}},$$

(45a)

$$\mathit{Q}.{\mathit{N}}_{\mathrm{C}}=0,\mathrm{on}{(\partial \mathcal{B})}_{\mathrm{C}},$$

(45b)

$${\mathit{Y}}_{\alpha }.{\mathit{N}}_{\mathrm{C}}=0,\mathrm{on}(\partial \mathcal{B}{)}_{\mathrm{C}},$$

(45c)

$$\hat{p}=0,\mathrm{on}{(\partial \mathcal{B})}_{\mathrm{left}}\mathrm{and}{(\partial \mathcal{B})}_{\mathrm{right}},$$

(45d)

$$\mathit{P}.{\mathit{N}}_A=0,\mathrm{on}{(\partial \mathcal{B})}_{\mathrm{left}}\mathrm{and}{(\partial \mathcal{B})}_{\mathrm{right}},$$

(45e)

where ${\mathit{N}}_{\mathrm{A}}=\pm {\mathit{E}}_Z$ is the unit normal to ${(\partial \mathcal{B})}_{\mathrm{left}}$ and ${(\partial \mathcal{B})}_{\mathrm{right}}$ , and ${\mathit{N}}_{\mathrm{C}}={\mathit{E}}_{\mathrm{R}}$ denotes the normal unit vector to ${(\partial \mathcal{B})}_{\mathrm{C}}$ . The boundary conditions specified in equations (45a)–(45e) assure that there is a constant concentration of nutrients at all times, the boundary ${(\partial \mathcal{B})}_{\mathrm{C}}$ is impermeable, there is no flux at the boundaries ${(\partial \mathcal{B})}_{\mathrm{left}}$ and ${(\partial \mathcal{B})}_{\mathrm{right}}$ , and the axial boundaries are free of pressure and traction, respectively, assuring that changes in the boundary are solely attributable to growth.

6.2. Axial mass flux

To comply with the assumptions made so far, we redefine the non-locality function so it depends only on z and $\tilde{z}$ . Therefore, by adapting the results obtained in Ramírez-Torres et al. [26] to our case, we set:

$$f_{\alpha (x-\tilde{x})}(x-\tilde{x}):=h_{\alpha (z,\tilde{z})}(z,\tilde{z})=\frac{L_{\mathrm{in}}^{\alpha (z,\tilde{z})-1}}{2\pi R_{\mathrm{in}}^2\mathrm{\Gamma }(1-\alpha (z,\tilde{z}))}{|z-\tilde{z}|}^{-\alpha (z,\tilde{z})},$$

(46)

where Γ(·) denotes Euler’s Gamma function and $\alpha (z,\tilde{z})\in ]0,1[$ . Consequently, the substitution of the last expression into equation (36) yields:

$${\left[{\mathit{y}}_{\alpha }(z,t)\right]}^z=-{\varrho }_f{\int }_{-\ell \left(t\right)}^{+\ell \left(t\right)}\frac{L_{\mathrm{in}}^{\alpha (z,\tilde{z})-1}}{2\mathrm{\Gamma }(1-\alpha (z,\tilde{z})){|z-\tilde{z}|}^{\alpha (z,\tilde{z})}}{\left[\mathit{d}(\tilde{z},t)\right]}^{\mathrm{zz}}\frac{\partial {\omega }_{\mathrm{cf}}(\tilde{z},t)}{\partial \tilde{z}}d\tilde{z},$$

(47)

where $\ell \left(t\right)>0$ denotes the changing axial length of the body. It should be noted that, in this work, we use a definition of mass flux vector, which is slightly different from the definition given in Ramírez-Torres et al. [26] (see equation (54) therein). Specifically, in Ramírez-Torres et al. [26], it was considered a fractional diffusivity tensor with anomalous physical units $L^{1+\alpha }{T}^{-1}$ whereas, in the present approach, we preferred to give the anomalous length units to the non-locality function with the aim of having a standard diffusivity tensor from the outset.

Remark 2 We observe that although the definition of the non-locality function presented in equation (46) bears similarities to the one proposed in Ramírez-Torres et al. [26], its motivation necessitates a different approach compared to the methodology followed in Ramírez-Torres et al. [26]. The primary challenge in addressing equation (47) lies in the difficulties arising in trying to express ${\left[{\mathit{y}}_{\alpha }(z,t)\right]}^z$ as a convolution integral. Nevertheless, we observe that by extending the domain to the real line, i.e., $\ell \to +\infty$ , and introducing the notation $y^z(\tilde{z},t):=-{\varrho }_f{\left[\mathit{d}(\tilde{z},t)\right]}^{\mathrm{zz}}\partial {\omega }_{\mathrm{cf}}(\tilde{z},t)/\partial \tilde{z}$ , we can reformulate equation (47) as:

$${\left[{\mathit{y}}_{\alpha }(z,t)\right]}^z={\int }_{-\infty }^{+\infty }\frac{L_{\mathrm{in}}^{\alpha (z,\tilde{z})-1}}{2\mathrm{\Gamma }(1-\alpha (z,\tilde{z})){|z-\tilde{z}|}^{\alpha (z,\tilde{z})}}\frac{\sin (\alpha (z,\tilde{z})\pi /2)}{\sin (\alpha (z,\tilde{z})\pi /2)}{y}^z(\tilde{z},t)d\tilde{z},$$

(48)

so that, by defining $I(z,\tilde{z}):=1/\mathrm{\Gamma }(1-\alpha (z,\tilde{z})){|z-\tilde{z}|}^{\alpha (z,\tilde{z})}\sin (\alpha (z,\tilde{z})\pi /2)$ and noticing that $\mathrm{\Gamma }\left(\alpha \right)\mathrm{\Gamma }(1-\alpha )=\pi /\sin \left(\alpha \pi \right)$ , we have that [107]:

$$I(z,\tilde{z})=\frac{\mathrm{\Gamma }\left(\alpha (z-\tilde{z})\right)}{\pi i^{\alpha (z,\tilde{z})}{(z-\tilde{z})}^{\alpha (z,\tilde{z})}}+\frac{\mathrm{\Gamma }\left(\alpha (z,\tilde{z})\right)}{\pi i^{\alpha (z,\tilde{z})}{(\tilde{z}-z)}^{\alpha (z,\tilde{z})}}=\frac{1}{\pi }{\int }_{-\infty }^{+\infty }{e^{i\xi (z-\tilde{z})}\left|\xi \right|}^{\alpha (z,\tilde{z})-1}\mathrm{d}\xi .$$

(49)

Thus, by substitution into equation (48), we can write:

$$\begin{array}{c}\multicolumn{1}{c}{{\left[{\mathit{y}}_a(z,t)\right]}^2=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{L}_{\mathrm{in}}^{\alpha {(z,\tilde{z})}^{-1}}\sin (\alpha (z,\tilde{z})\pi /2)\left\{{\int }_{-\infty }^{+\infty }{e}^{i\xi (z-\tilde{z})}{\left|\xi \right|}^{\alpha (z,\tilde{z})-1}\mathrm{d}\xi \right\}{y}^z(\tilde{z},t)d\tilde{z}}\\ \multicolumn{1}{c}{=\frac{1}{2\pi }{\int }_{-\infty }^{+\infty }{e}^{1\xi z}\left\{{\int }_{-\infty }^{+\infty }\left[L_{\mathrm{in}}^{\alpha (z,\tilde{z})-1}\sin (\alpha (z,\tilde{z})\pi /2)\right]{\left|\xi \right|}^{\alpha (z,\tilde{z})-1}{e}^{-i\xi \tilde{z}}{y}^z(\tilde{z},t)\mathrm{d}\tilde{z}\right\}\mathrm{d}\xi .}\end{array}$$

(50)

Hence, when $\alpha$ approaches either one or zero, and z and $\tilde{z}$ are fixed, it can be inferred that:

$$\underset{\alpha \to 1^{-}}{lim}{\left[{\mathit{y}}_{\alpha }(z,t)\right]}^z={\mathcal{F}}^{-1}\{\mathcal{F}\{y^z\}(\xi ,t)\}(z,t)=y^z(z,t)=-{\varrho }_{\mathrm{f}}{\left[\mathit{d}(z,t)\right]}^{\mathrm{zz}}\frac{\partial {\omega }_{\mathrm{cf}}}{{\partial }_z}(z,t),$$

(51a)

$$\underset{\alpha \to 1^{-}}{lim}{\left[{\mathit{y}}_{\alpha }(z,t)\right]}^z=0,$$

(51b)

where $\mathcal{F}\{·\}$ denotes the Fourier transform and ${\mathcal{F}}^{-1}\{·\}$ denotes its inverse. Thus, although a different approach was used in Ramírez-Torres et al. [26] to formulate ${\mathit{y}}_{\alpha }$ , the findings in equations (51a) and (51b) still maintain the original non-local diffusion patterns proposed in Ramírez-Torres et al. [26]. These, in particular, enable the examination of the connection between the function $\alpha$ and the dynamics of nutrient diffusion. When $\alpha$ approaches one, it corresponds to the conventional Fick’s law. Conversely, when $\alpha$ approaches zero, it signifies that the diffusion of nutrients is impeded, reflecting hindered transport.

The restrictions imposed on the motion, together with the definitions (38a) and (46), imply that the material non-locality function can be rephrased as:

$${\mathfrak{h}}_{\stackrel{⌣}{a}(z,\tilde{z},t)}(z,\tilde{z},t)=\frac{[1+u\prime (z,t)]L_{\mathrm{in}}^{\stackrel{⌣}{a}(z,\tilde{z},t)}}{2\pi R_{\mathrm{in}}^2\mathrm{\Gamma }(1-\stackrel{⌣}{\alpha }(Z,\tilde{Z},t))}\frac{1}{{|Z+u(Z,t)-\tilde{Z}-u(\tilde{Z},t)|}^{\stackrel{⌣}{a}(z,\tilde{z},t)}},$$

(52)

where the notation $u\prime (Z,t)$ has been used to denote the partial derivative of u with respect to Z and:

$$\stackrel{⌣}{\alpha }(Z,\tilde{Z},t):=\alpha ({\chi }^Z(z,t),{\chi }^Z(\tilde{z},t))=\alpha (Z+u(Z,t),\tilde{Z}+u(\tilde{Z},t)).$$

(53)

Furthermore, the axial component of the material diffusivity tensor given in equation (38b) is:

$${\left[\mathfrak{D}(Z,\tilde{Z},t)\right]}^{\mathrm{ZZ}}=\frac{1+u\prime (\tilde{Z},t)-{\gamma }^3(\tilde{Z},t){\mathrm{\Phi }}_{\mathrm{s}v}}{[1+u\prime (Z,t)][1+u\prime (\tilde{Z},t)]}{d}_{\mathrm{R}}.$$

(54)

Thus, considering equations (52) and (54), the axial component of the material flux vector ${\mathit{Y}}_{\alpha }$ specified in equation (37) can be written in the form:

$${\left[{\mathit{Y}}_{\alpha }(Z,t)\right]}^Z=-{\varrho }_f{\int }_{-L_{\mathrm{in}}}^{+L_{\mathrm{in}}}\frac{[1+u\prime (\tilde{Z},t)-{\gamma }^3(\tilde{Z},t){\mathrm{\Phi }}_{\mathrm{s}v}]{[1+u\prime (\tilde{Z},t)]}^{-1}{L}_{\mathrm{in}}^{\stackrel{⌣}{a}(Z,\tilde{Z},t)-1}{d}_{\mathrm{R}}{\hat{\omega }\prime }_{\mathrm{cf}}(\tilde{Z},t)}{2\mathrm{\Gamma }(1-\stackrel{⌣}{a}(Z,\tilde{Z},t)){|Z+u(Z,t)-\tilde{Z}-u(\tilde{Z},t)|}^{\stackrel{⌣}{a}(Z,\tilde{Z},t)}}\mathrm{d}\tilde{Z}.$$

(55)

7. Numerical results

In this section, we present the numerical results obtained through the implementation in COMSOL Multiphysics® of the model equations in their weak formulation. These results will serve as a critical tool in our quest to unravel the complexities of our model and shed light on its implications. In particular, the parameters used in our simulations are reported in Table 1.

Table 1.

List of parameters used in the numerical simulations.

Parameter	Unit	Value	Parameter	Unit	Value
$L_{\mathrm{in}}$	cm	0.500	${\varrho }_{\mathrm{s}}$	kg/m3	1000
$R_{\mathrm{in}}$	cm	1.000 x10—2	${\varrho }_{\mathrm{f}}$	kg/m3	1000
$\lambda$	Pa	1.333x104	$d_{\mathrm{R}}$	m2/s	3.200 x 10—9
$\mu$	Pa	1.999 x 104	${\zeta }_{\mathrm{fp}}$	kg/(m3 s)	1.343 x 10—3
$k_{\mathrm{R}}$	m2/(Pa s)	4.875 x 10—13	${\zeta }_{\mathrm{nf}}$	kg/(m3 s)	1.150 x 10—5
$m_0$	—	0.0848	${\zeta }_{\mathrm{cp}}$	kg/(m3 s)	3.000 x 10—4
$m_1$	—	4.638	${\zeta }_{\mathrm{pn}}$	kg/(m3 s)	1.500 x 10—3
${\delta }_1$	—	7.138 x 10—1	${\omega }_{\mathrm{cr}}$	—	1.000 x 10—3
${\delta }_{\mathit{2}}$	Pa	1.541 x 103	${\omega }_{\mathrm{env}}$	—	7.000 x 10—3
${\mathrm{\Phi }}_{\mathrm{s}v}$	—	0.8	${\omega }_0$	—	1.000 x 10—2

The values were taken from different sources [9,18,24,26,87,97,108,109].

We notice that, in addition to the parameters listed in Table 1, we also need to consider the Helmholtz free energy of the fluid phase. Following Penta et al. [110], we consider the mass fraction of nutrients to be low, so that ${\hat{\psi }}_{\mathrm{f}}$ can be approximated by [111, 110]:

$${\stackrel{⌣}{\psi }}_f=C_f[f\left(T\right)-T\log (T/T_0)],$$

(56)

where C_f is the specific heat capacity set to be equal to the specific heat of water 4184 J/(kg K), f (T) is a function of temperature, assumed here to be linear in T with slope b > 0, and T₀ denotes a reference value for the temperature. Since we are dealing with isothermal processes, T is constant and, for simplicity, we assume that T=T₀=310.15 K. Moreover, we set b=10^–6.

7.1. Case I: constant-order non-local flux

As a first step in our investigations, we begin by assuming a constant fractional parameter and explore how the new growth law (refer to equation (44c)) impacts the dynamics of an avascular tumour subjected to a non-Fickian diffusion law. Essentially, our analysis focuses on understanding the influence of $\mathcal{H}$ on the tumour’s evolution. To achieve this, we assign ${\kappa }_{\mathrm{fp}}$ and ${\kappa }_{\mathrm{nf}}$ the same value of $a_v=2\times {10}^{11}$ Pa s, so that we obtain the same growth law as considered in Ramírez-Torres [26] when $\mathcal{H}$ is absent.

Figure 1 presents the spatial distributions of the axial displacement u and the growth law $R_{\gamma }$ at a fixed time of 20 days. The comparison in Figure 1 focuses on the model’s outcomes by switching on and off the contribution of the configurational mechanical stress $\mathcal{H}$ for two different values of the fractional parameter $\alpha$ . Specifically, the fractional parameter is set equal to $\alpha =0.1$ , indicating hindered diffusion, and to $\alpha =0.9$ , a scenario of nearly Fickian diffusion (see Remark 2). By examining the left panel of Figure 1 for $\alpha =0.9$ , we observe that the displacement is higher in the model that neglects $\mathcal{H}$ with respect to the model that accounts for this term. This behaviour can be explained by referring to the right panel of Figure 1 as the sign of $R_{\gamma }$ determines if there is mass accretion $(R_{\gamma }>0)$ or mass resorption $(R_{\gamma }<0)$ and its value is directly related to the amount of displacement the tumour experiences. In fact, in the case without $\mathcal{H}$ , $R_{\gamma }$ becomes negative when the nutrients’ mass fraction goes below the critical value, namely, ${\hat{\omega }}_{\mathrm{cf}}<{\omega }_{\mathrm{cr}}$ (as shown in Figure 2). However, the sign of $R_{\gamma }$ can also change when ${\hat{\omega }}_{\mathrm{cf}}>{\omega }_{\mathrm{cr}}$ due to the contribution of $\mathcal{H}$ being larger than $\mathrm{\Upsilon }:=3[{\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}}]/{\alpha }_v$ (refer to equation (44c)). This conclusion is supported by both the right panel of Figure 1 and the left panel of Figure 2. In the latter, we present the spatial distribution of $3\mathcal{H}/\left(a_v{J}_{\gamma }\right)$ and $\mathrm{\Upsilon }$ for $\alpha =0.9$ . Indeed, we notice that when the configurational stress is accounted for, the growth law becomes negative at around Z = 0.24 cm, while the nutrients’ mass fraction falls below the critical value (represented with the dashed horizontal line in the right panel of Figure 2) at around Z = 0.08 cm. Therefore, in our model, mass resorption can still occur in the body even if the nutrients’ mass fraction does not fall below the critical value. This behaviour showcases the influence of the configurational mechanical stress on the body’s dynamics and its direct impact on growth, as per our adopted approach. While the growth law without $\mathcal{H}$ relies heavily on the presence and changes in nutrient availability, the current model incorporates the Eshelby stress in the growth law, which, for $\alpha =0.9$ , prevents the displacement from increasing. In contrast, if $\alpha =0.1$ , we observe that the numerical results for the case that includes the contribution of $\mathcal{H}$ (represented by the dotted line with square markers) exhibit an increase in the axial displacement that is less concentrated at the axial boundary compared to the model that neglects $\mathcal{H}$ (represented by the solid line with plus sign markers). By referring to the right panel of Figure 1, we notice that while $R_{\gamma }<0$ in almost the whole spatial domain in the absence of $\mathcal{H}$ and $\alpha =0.1$ , it is positive in large part of the spatial domain when $\mathcal{H}$ is not neglected and reflects the contribution to growth given by the configurational mechanical stress by aiding the displacement to increase and not be confined to the axial boundary.

Figure 1.

(Left panel): spatial distribution of the axial displacement u. (Right panel): spatial distribution of the growth law $R_{\gamma }$ . The model’s outcomes are compared by toggling the contribution of $\mathcal{H}$ . All curves shown represent the results at a fixed time of 20 days.

Figure 2.

(Left panel): spatial distribution of $3\mathcal{H}/\left(a_v{J}_{\gamma }\right)$ and $\mathrm{\Upsilon }:=3({\kappa }_{\mathrm{fp}}{\mathrm{\Upsilon }}_{\mathrm{fp}}-{\kappa }_{\mathrm{nf}}{\mathrm{\Upsilon }}_{\mathrm{nf}})/{\alpha }_v$ for $\alpha =0.9$ . (Right panel): spatial distribution of the nutrients’ mass concentration ${\hat{\omega }}_{\mathrm{cf}}$ . All curves shown represent the results at a fixed time of 20 days. The horizontal dashed line represents the critical value ${\omega }_{\mathrm{cr}}$ .

To continue with our analysis, in Figure 3, we present the spatial profiles of the pressure and the mass fraction of the proliferating cells. Focusing on the left panel of Figure 3, we observe that necrotic cells appear when ${\hat{\omega }}_{\mathrm{cf}}<{\omega }_{\mathrm{cr}}$ , which is consistent with the spatial distribution of nutrient mass fraction shown in Figure 2. Particularly, in examining the case in which $\mathcal{H}$ is active and $\alpha =0.9$ , the model foresees a low mass fraction of necrotic cells (although not visible to the unaided eye) compared to the scenario without $\mathcal{H}$ , where a larger necrotic cell’s mass fraction is appreciated. However, when $\mathcal{H}$ is present, we note that the displacement is lower than in the case without $\mathcal{H}$ . This behaviour emphasises how the properties of ${\mathit{F}}_{\gamma }$ play a major role in the equation that governs cell proliferation through the term r_ph. Still, the appearance of the necrotic region is ascribable to ${\hat{\omega }}_{\mathrm{cf}}$ being less than ${\omega }_{\mathrm{cr}}$ .

Figure 3.

(Left panel): spatial distribution of the proliferating cells’s mass fraction ${\hat{\omega }}_{\mathrm{ps}}$ . (Right panel): spatial distribution of the pressure $\hat{p}$ . The model’s outcomes are compared by toggling the contribution of $\mathcal{H}$ . All curves shown represent the results at a fixed time of 20 days.

Now, focusing on the right panel of Figure 3, we observe that in the absence of $\mathcal{H}$ and for $\alpha$ = 0.1, the pressure is positive. This phenomenon can be attributed to hindered diffusion resulting from selecting $\alpha$ = 0.1, which leads to the existence of a large region of necrotic cells at the centre of the tumour. In agreement with the biological foundations of nutrient diffusion and necrosis in a tumour [26,112], the dissolution of necrotic cells in water generates an outward flux, opposing the direction of fluid flow, which makes p to be positive. We notice that the red solid curve undergoes a change in monotony when a = 0.9, and this phenomenon can be attributed to the increasing prominence of the mass fraction of necrotic cells. However, when considering the presence of H and a = 0.1, a significant portion of necrotic cells no longer generates a positive pressure. Furthermore, if a = 0.9, for which the concentration of necrotic cells is low, the pressure becomes positive. In such a case, the change of sign of the pressure (as it is negative for initial times) appears to be strongly associated with the configurational mechanical stress rather than the evolution of the necrotic cells’ mass fraction as the coupling of the dynamic equations with the latter is weaker.

In the left panel of Figure 4, we illustrate the temporal evolution of the tumour’s volume across a 100-day span. We remark that, our model, built upon the specific choice of the scalar permeability (refer to equation (21)), necessitates that $J/J_{\gamma }>{\mathrm{\Phi }}_{\mathrm{s}v}$ for all $X\in \mathcal{B}$ and $t\in \mathcal{I}$ . Figure 4 demonstrates that in the absence of the configurational mechanical stress in the growth law, this criterion is breached at the axial boundary before the specimen’s volume stabilises (see right panel of Figure 4). Specifically, for $\alpha =0.1$ , the ratio $J/J_{\gamma }$ reaches ${\mathrm{\Phi }}_{\mathrm{s}v}$ when evaluated at $Z=L_{\mathrm{in}}$ at around approximately 55 days. Furthermore, for $\alpha =0.9$ , it occurs after about 85 days. As a result, our simulations, constrained by the specific parameter selections and underlying constitutive assumptions, are unable to extend beyond these time frames. Notably, both scenarios without $\mathcal{H}$ exhibit a phase of resorption and, in particular, for $\alpha =0.1$ , the tumour reverts to its initial volume shown with the horizontal dashed line in the left panel of Figure 4. Conversely, when $\mathcal{H}$ is accounted for, the specimen’s volume stabilises while maintaining $J/J_{\gamma }$ greater than the volumetric fraction of the solid phase in the natural state, ${\mathrm{\Phi }}_{\mathrm{s}v}$ , for all $X\in \mathcal{B}$ and $t\in \mathcal{I}$ . These considerations highlight the pivotal role of the configurational stress in our model and its influence on the evolution of the tumour.

Figure 4.

(Left panel): time evolution of the tumour’s volume. (Right panel): spatial distribution of the ratio $J/J_{\gamma }$ at different instants of time.

Before going further, we study the influence of a_v in the evolution of the body specimen. Specifically, in Figure 5, distinct growth patterns can be observed for different values of a_v. It is noteworthy that for $\alpha$ = 0.9, as a_v increases, the displacement distribution tends to resemble the case where the contribution of H in the growth law is absent. In addition, when $\alpha$ = 0.1, decreasing a_v leads to a less pronounced impact of nutrients at the axial boundary. This is evident from the nearly identical curves for $a_v=1\times {10}^8$ Pa s and $a_v=1\times {10}^{10}$ Pa s, indicating that the coupling with the configurational stress is stronger compared to the coupling with the evolution of nutrients.

Figure 5.

(Left panel): spatial distribution of the axial displacement for $\alpha =0.9$ . (Right panel): spatial distribution of the axial displacement for $\alpha =0.1$ . The model’s outcomes are compared by selecting different values of a_v. All curves shown represent the results at a fixed time of 20 days.

7.2. Case II: variable-order non-local flux

In this second scenario, we explore how a variable fractional-order affects the evolution of the main parameters in our model. Before delving deeper into the analysis, we note that although there is a lack of experimental data, we opt for a functional form that holds a physical interpretation. Specifically, we prescribe $\alpha$ in a way that it acts as an overall descriptor of the complex micro-environment found in tumours by reflecting, for instance, a poorly organised vasculature leading to regions of limited diffusion [69]. To facilitate our study, we consider that $\alpha (z,\tilde{z})={\alpha }_1\left(\tilde{z}\right)$ and set:

$${\alpha }_1\left(\tilde{z}\right):=A-B\exp (-\frac{1}{C^2}{(\tilde{z}-L_{\mathrm{in}}/2)}^2),$$

(57)

where A, B, and C are positive constants, with A and B being dimensionless and C having physical units of length. As depicted in Figure 6, at the initial time and for the given values of A, B, and C, the expression in equation (57) enables the description of a non-diffusive regime in the vicinity of $\pm L_{\mathrm{in}}/2$ . Therefore, our focus lies in understanding how this particular selection influences the tumour dynamics.

Figure 6.

Profile of the variable-order parameter at t = 0 s and with A = 0.9, B = 0.85, and $C=2\times {10}^{-2}$ cm.

Figure 7 illustrates the evolution of the spatial profiles of the axial displacement and the mass concentration of nutrients. We notice that the variable and non-local way in which the nutrients diffuse into the tissue affects the manner in which the tumour grows. Specifically, we observe a notable inhomogeneous spatial distribution in the axial displacement, which corresponds to a significant decrease in the nutrients’ mass concentration near $Z=L_{\mathrm{in}}/2$ as prescribed by equation (57). In particular, even though the variable-order ${\alpha }_1$ returns to 0.9 in the interval $[0,L_{\mathrm{in}}/2]$ , which corresponds to an unhindered diffusion region, the nutrient availability remains low. This outcome can be attributed to our decision to prescribe $\alpha$ as a function of $\tilde{z}$ , giving importance to a situation where $\alpha$ retains a stronger “spatial memory” of its previous states (see Remark 1). We emphasise that, in contrast, the case in which $\alpha$ depends solely on z would result in disregarding information from remote points and causing the mass flux to respond relatively quickly to changes. Such an approach would require facing challenges both from numerical and physical points of view, similar to those observed in previous studies (refer to Sun et al. [113] for an example).

Figure 7.

(Left panel): spatial distribution of the axial displacement u. (Right panel): spatial distribution of the nutrients’ mass concentration ${\hat{\omega }}_{\mathrm{cf}}$ . The curves correspond to the fixed times of t = 1 day and t = 20 days.

To continue our analysis, we present the numerical results for the mass fraction of proliferating cells and the growth parameter. The left panel of Figure 8 shows a significant decline in proliferating cells where ${\alpha }_1$ reaches its minimum value. This minimum value corresponds to the spatial location with the greatest inhibition to diffusion, creating a region of space of limited nutrient availability. The right panel of Figure 8 depicts the spatial profile of the growth parameter γ, which exhibits a similar sudden “jump” at $L_{\mathrm{in}}/2$ . This jump increases over time in correspondence with the axial displacement evolution. We also notice that the part of the domain in which the necrotic cells appear coincides with the one in which the nutrients fall below the critical value. Furthermore, when examining Figure 9, we can observe a decrease in the pressure of the interstitial fluid over time from the free boundary towards the centre of the tumour which, in the case reported here, does not become positive. By making reference to our discussions in the previous section, this can be attributed to two factors: first, the region of necrotic cells is not sufficiently large to generate an outward flux, and second, the influence of $\mathrm{\Upsilon }$ on the pressure is dominant, surpassing the contribution of $3\mathcal{H}/\left(a_v{J}_{\gamma }\right)$ , as evidenced by the positive values of $R_{\gamma }$ in the right panel of Figure 9.

Figure 8.

(Left panel): spatial distribution of the proliferating cells’ mass fraction ${\hat{\omega }}_{\mathrm{ps}}$ . (Right panel): spatial distribution of the proliferating cells’ mass fraction. The curves correspond to the fixed times of t = 1 day and t = 20 days.

Figure 9.

(Left panel): spatial distribution of the pressure. (Right panel): spatial distribution of the growth law. The curves correspond to the fixed times of t = 1 day and t = 20 days.

An important observation arises from the spatially inhomogeneous growth giving rise to residual stresses (see, for example, [13, 15, 39, 114–117] and references therein), i.e., the stress in the material is different from zero even though there are no applied loads. As reported in Figures 1 and 7, the inhomogeneous nature of the displacement suggests the existence of regions where the tumour grows faster due to the non-uniform distribution of nutrients. Particularly, referring to Figure 7, this phenomenon manifests in the region $[L_{\mathrm{in}}/2,L_{\mathrm{in}}]$ (and, by symmetry, in $[-L_{\mathrm{in}},-L_{\mathrm{in}}/2]$ ) where both the proliferating cells’ mass fraction, ${\hat{\omega }}_{\mathrm{ps}}$ , and the growth parameter, γ, exhibit higher values (see Figure 8). As illustrated on the left panel of Figure 10, the circumferential stress (and, due to the inherent symmetry, the radial stress as well) is more compressive in the spatial region where the specimen experiences increased growth, which is indicative of the proliferation activity. On the contrary, in the inner region (i.e., $[0,L_{\mathrm{in}}/2[$ ), where the tumour exhibits reduced growth, the circumferential stress is less compressive. Although experimental and theoretical studies are mainly focused on the growth of tumour spheroids [9, 13, 108, 116, our theoretical findings align with evidence supporting the existence of residual stresses induced by growth and of stresses compressively increasing within the tumour region. We remark that, in the context of this work, we are not considering the restriction by a surrounding tissue.

Figure 10.

Spatial distribution of the circumferential (left panel) and axial (right panel) stresses. The curves correspond to the fixed times of t = 1 day and t = 20 days.

If we direct our attention to the right panel in Figure 10, the axial stress vanishes at $Z=L_{\mathrm{in}}$ and, at initial times, is compressive. Particularly, in the neighbourhood of the point $Z=L_{\mathrm{in}}/2$ , where diffusion is hindered (see equation (57) and Figure 7), the axial stress exhibits a non-monotonic trend and, as time progresses, a jump-like transition can be observed. To gain insight into this pattern, we notice that although the decline in nutrient availability becomes less pronounced over time, it eventually falls below the critical threshold ${\omega }_{\mathrm{cr}}$ . As reported in Figure 8, this event effectively creates a region of actively proliferating cells and another one with presence of necrotic cells, which underscores the sharp jump of the axial stress at $Z=L_{\mathrm{in}}/2$ and, consequently, the critical role of the non-uniform distribution of nutrients in shaping the tumour’s dynamics. We further notice that at points near $Z=L_{\mathrm{in}}/2$ from the right, the stress is in tension, whereas it is in compression as we approach $Z=L_{\mathrm{in}}/2$ from the left. This trend bears a resemblance to patterns documented in studies examining tumour growth in contact with host tissues (as in, for example, [9, 108, 118]). In our specific scenario, we attribute this behaviour to both the sharp region of restricted diffusion and the abrupt alteration in the tumour’s micro-environment (refer to Figure 8).

8. Conclusion

This study introduces a novel approach to model the growth of an avascular tumour by incorporating the effects of nutrient transport through a non-local constitutive relationship. Notably, we extend the methodology in Ramírez-Torres et al. [26] by proposing a non-local law for the mass flux vector that introduces a non-uniform parameter to include relevant information about the complex tumour micro-environment, such as regions with hindered diffusion. Our framework has similarities with the notion of variable-order fractional operators, which has been acknowledged for their potential in simulating interdisciplinary processes [58–61]. It is important to note that both the non-locality of the constitutive law and the variable-order parameter evolve with the growth-induced inelastic distortions that occur during the evolution of the system and its visible deformation. Hence, by allowing the variable-order parameter to vary, our approach offers a potential framework for capturing space-dependent microstructural phenomena with the variable-order parameter $\alpha$ serving as the link transferring the information across the scales.

In our simulations, we explored the impact of non-Fickian diffusion on tumour progression in a simplified yet significant scenario involving the axial growth of a cylindrical specimen. Within this framework, owing to imposed symmetries while still recognising its three-dimensional nature, we consider a non-locality function, $f_{\alpha }$ , which captures the dependence of the non-local mass flux vector on both elastic deformation and inelastic distortions (refer to equation (55)). Our choice of the non-locality function is akin to the one presented in Ramírez-Torres et al. [26], with the distinction that our approach incorporates a non-uniform variable-order fractional operator (refer to equation (46)). Importantly, this modification maintains a key characteristic of the definition in Ramírez-Torres et al. [26]. As discussed in Remark 2, when $\alpha$ tends to one, the mass flux vector aligns with conventional Fick’s law. Conversely, as $\alpha$ approaches zero, the diffusion of nutrients becomes hindered. In this respect, our model introduces a notable improvement to enhance realism by considering the effects of fractional parameter inhomogeneity. We hypothesise that the fractional parameter $\alpha$ characterises, in a homogenised sense, situations where the tumour is composed of areas with varying diffusion properties. Specifically, we examine the transition to a region with low diffusion in an attempt to establish a correlation between the fractional parameter and the micro-environment of the tumour tissue, as depicted in Figure 6.

In contrast to the constant-order non-local flux model, the variable-order non-local model may acknowledge the complexity of the medium by accounting for spatial phenomena arising in a heterogeneous and complex environment. This complexity, as per our proposed model, is characterised by the presence of regions with low diffusion within the tumour. Our numerical simulations revealed that the evolution of the main parameters in the model is significantly influenced by this selection and, notably, we observed more complex behaviours in displacement, growth, nutrients’ mass fraction, the distribution of proliferating cells, and stresses. Overall, these findings provide valuable insights into the complex dynamics and highlight the potential of applying this modelling framework to study even more intricate scenarios in the future.

Another crucial aspect of our work, distinguishing it from Ramírez-Torres et al. [26], is the establishment of a new growth law following the approach presented in Grillo and Di Stefano [70, 71]. Within this setting, the evolution law for the inelastic distortions is determined without relying on initial phenomenological assumptions. As a consequence, we introduced non-standard source/sink terms in the mass balance equations for proliferating and necrotic cells accounting for the effects of the configurational mechanical stress. These terms are determined once the equation for growth-induced inelastic distortions has been established through the study of the dissipation inequality. To account for the implications of embracing this approach, we investigated a benchmark scenario under the assumption of a uniform fractional parameter, considering two possibilities: the presence or absence of the term $\mathcal{H}$ in the growth law. In particular, incorporating $\mathcal{H}$ has a significant impact on the spatial distribution and temporal evolution of the model’s unknowns with respect to when it is ignored. The inclusion of the configurational mechanical stress in the growth law predicts that the accumulation and removal of mass cannot be solely attributed to nutrient availability thereby H reflects a positive or negative contribution to growth.

Our work holds potential for future enhancements and broader applicability. One area pertains to tailoring the mass flux vector of nutrients to the dimensions and symmetries of the problem. Moreover, further exploration is needed to understand the connection of the variable-order parameter to the tumour micro-environment or other physical phenomena that may influence it. Therefore, it is crucial to develop constitutive relationships in collaboration with chemists and biologists for future research. Establishing such relationships could yield valuable insights in the analysis of drug transport and treatment within tumour tissues. As recognised in different studies (see, for example, [119]), the tumour’s complex internal structure plays a crucial role in controlling how drugs disperse through soft tissues. Among several factors, the chaotic and aberrant nature of tumour blood vessels, characterised by tortuosity and irregular connections disrupting blood flow, hampers the performance of anti-tumour drugs [120]. Thus, accounting for the spatial aspects of diffusion through the tumour’s micro-environment and their influence on drug delivery in mathematical models is highly desirable to enhance therapeutic efficacy. The framework proposed in this work has the potential to address the spatial complexities of drug transport within tumour tissues. Indeed, in our model, the introduction of an integro-differential operator of variable-order offers a mathematical tool to encode, in an effective way, spatial irregularities that change with the dynamics of the tissue and their impact on distant points. Part of our current efforts aim to pinpoint areas where drug diffusion could potentially be insufficient by deciphering how the variability in $\alpha$ may be linked to inconsistent and less efficient drug delivery. In exploring the complex interplay between vascular abnormalities and drug diffusion, our model might suggest a pathway toward a more detailed understanding of drugs’ non-local diffusion mechanisms due to the heterogeneous tumour’s micro-environment. This knowledge might contribute to developing improved drug delivery strategies, which could, in turn, enhance treatment outcomes.

After witnessing significant advancements in fractional calculus across various scientific disciplines over the past decade, the incorporation of such a mathematical tool into the field of biological sciences could yield remarkable contributions. While the use of fractional calculus in real-world phenomena is still in its early stages, it has proven to be an exceptional mathematical tool for comprehending extraordinary events. Our model endeavours to capture the intricate complexities observed in real tumours, striving to replicate the complex tumour growth patterns and the distribution of nutrients within the tumour micro-environment with a level of realism. In doing so, it offers valuable insights into the dynamics of avascular tumour growth and emphasises the significance of considering the non-local effects of the tumour micro-environment in tumour growth models.