A Genuinely Hybrid, Multiscale 3D Cancer Invasion and Metastasis Modelling Framework

We introduce in this paper substantial enhancements to a previously proposed hybrid multiscale cancer invasion modelling framework to better reflect the biological reality and dynamics of cancer. These model updates contribute to a more accurate representation of cancer dynamics, they provide deeper insights and enhance our predictive capabilities. Key updates include the integration of porous medium-like diffusion for the evolution of Epithelial-like Cancer Cells and other essential cellular constituents of the system, more realistic modelling of Epithelial–Mesenchymal Transition and Mesenchymal–Epithelial Transition models with the inclusion of Transforming Growth Factor beta within the tumour microenvironment, and the introduction of Compound Poisson Process in the Stochastic Differential Equations that describe the migration behaviour of the Mesenchymal-like Cancer Cells. Another innovative feature of the model is its extension into a multi-organ metastatic framework. This framework connects various organs through a circulatory network, enabling the study of how cancer cells spread to secondary sites.


Introduction
Mathematical modelling has been a powerful tool in cancer research, aiding in the understanding of the intricate biological processes, in predicting cancer growth and metastasis, and in the treatment of the disease.However, the very dynamic nature of cancer, characterized by tumour heterogeneity, complex intercellular interactions, and evolutionary dynamics, presents challenges for traditional mathematical models.A particularly challenging aspect of cancer growth lies in the understanding of the intricate interplay between epithelial and mesenchymal phenotypes, which can greatly influence tumour progression and response to therapy.
In this context, it was previously proposed in Sfakianakis et al. (2020), Franssen et al. (2021), a hybrid mathematical modelling framework capable of capturing the dynamic transitions between these phenotypes and shedding light on their underlying mechanisms.By combining deterministic and stochastic mathematical modelling frameworks, this hybrid approach is able to account for the heterogeneous nature of cancer cells, the complex macroscopic dynamics, and the stochasticity inherent in cellular migration.This integrated approach enables more comprehensive modelling of the Epithelial-Mesenchymal Transition (EMT), that enables the individual cell invasion, and of the reverse procedure of Mesenchymal-Epithelial Transition (MET) dynamics in cancer cells, cf.Roche (2018), offering insights into critical processes such as local tissue invasion, tumour island formation, and cancer metastasis.This previously proposed model and work have provided with valuable insights into the understanding of the relevant biological processes and the capabilities of the employed mathematics.
There are relatively few studies in the literature that share a similar modelling philosophy.Notable examples include Colombi et al. (2015a), Colombi et al. (2015b), where a measure theoretic approach was aplpied to explore cell differentiation and aggregation; Colombi et al. (2017), where a biological and mathematical "switch" between concentrated cell-particle descriptions and distributed mass approaches was invesitgated; and Capasso and Morale (2010), where a "doubly stochastic" system of interacting cell-particles was analyzed, demonstrating that in the limit of a large number of cell-particles, deterministic Partial Differential Equations (PDEs) emerge.Moreover, Te Boekhorst et al. (2016) discusses how cell migration, influenced by mechanical and chemical interactions with the extracellular environment, is modulated by cellular and tissue determinants, affecting decision-making and migration outcomes in processes like morphogenesis, repair, immune surveillance, and cancer metastasis.Further contributions include Hiremath and Surulescu (2016), Hiremath et al. (2018) , Colombi and Scianna (2017), who examined coupled systems of PDEs and stochastic differential equations (SDEs) to describe biological processes across population, cellular, and sub-cellular levels.Additionally, Cañizo et al. (2015); Carrillo et al. (2018) focused on continuum models for interacting cell-particles, establishing the existence of global minimizers for these systems.
In the current paper, we propose several model updates that better align with the accumulated understanding, incorporate new methodologies, and address the limitations of the original modelling framework, ultimately paving the way for even deeper insights, improved cancer predictions, and broader applicability in cancer modelling.The model extensions that we propose, incorporate among others porous medium-like diffusion for the time evolution of the Epithelial-like Cancer Cells (ECCs) and the various other living cell components of the system.Further model extensions account for more biologically realistic EMT and MET modelling by accounting for the role of Cancer-Associated Fibroblast (CAF) cells and Transforming Growth Factor beta (TGF-β) in the tumour microenvironment.
CAFs are key players in the tumor microenvironment, affecting cancer progression, metastasis, and even therapy response.Their diverse roles, from the remodelling of the ECM to the modulation of growth factor signaling, facilitate tumor growth and metastasis, emphasizing the complexity of their interactions with cancer cells Weinberg 2014; Erdogan and Webb 2017.CAFs drive EMT in cancer cells through paracrine TGF-β signaling, a process crucial for metastasis, Yu et al. (2014), Xing (2010), Darby et al. (2014) further illustrate CAFs' contributions to the supportive tumor environment and parallel their functions in wound healing, highlighting the pathological activation of fibroblasts, i.e. of CAFS, in cancer.This highlights the significance of CAFs in the Biology of cancer and the potential of targeting their signaling pathways for therapeutic intervention.
TGF-β plays a dual role in cancer progression, acting initially as a tumor suppressor by inhibiting cell proliferation and inducing apoptosis, and subsequently promoting metastasis through EMT.This transition, crucial for both embryonic development and cancer metastasis, involves epithelial cells acquiring mesenchymal traits, enhancing their invasiveness and resistance to apoptosis, see e.g.Bierie and Moses (2006), Kalluri and Weinberg (June 2009).In Oft et al. (1998), TGF-β's essential role in tumor cell invasiveness was established, with further studies, Xu et al. (2009), Katsuno et al. (2013), detailing its mechanisms.Later, in Lee and Massagué (2022), further highlighted TGF-β's involvement in tissue remodeling and fibrosis, extending its impact beyond cancer to the tumor microenvironment.This body of work underscores TGFβ's complex role in caner progression and EMT in particular.
In addition, the model features updated Stochastic Differential Equations (SDEs) for the migration of Mesenchymal-like Cancer Cells (MCCs).Most importantly though, we present here for the first time in our modelling approach, the progress we have made using our hybrid model on a multiple organ metastatic framework.
Mathematical modelling of cancer invasion goes back almost thirty years to the work of Gatenby and Gawlinski (1996) and Perumpanani and colleagues (1996).These models were systems of nonlinear reaction-diffusion-taxis equations and studied travelling wave solutions modelling the invasive cancer cells cf.Marchant et al. (2000), Gatenby et al. (2006).Subsequent models considered hybrid discrete-continuum approaches (Anderson et al. 2000), various discrete or individual-based approaches (Turner and Sherratt 2002;Smallbone et al. 2005), nonlocal modelling focussing on cell-cell and cell-matrix adhesion (Domschke et al. 2014) and different multiscale approaches (Chaplain and Lolas 2005;Shuttleworth and Trucu 2019a, b, c;Franssen et al. 2019).
For a more comprehensive overview of the previous cancer modelling literature, the reader is referred to the recent review article of Sfakianakis and Chaplain (2021).
In identifying metastases, the integration of mathematical models is of significant importance as it can potentially offer a predictive approach that can complement existing technologies.Mathematical models could enhance the precision and efficiency of detecting micrometastatic sites (Banyard and Bielenberg 2015;Zavyalova et al. 2019;Chiang et al. 2016;Fares et al. 2020).It is also expected that these advanced predictive capabilities will be crucial for tailoring surgical and therapeutic strategies (Smith et al. 2019), potentially transforming the landscape of preoperative and intraoperative intervention.
The rest of the paper is structured as follows: in Sect. 2 we briefly describe the previous mathematical model, in Sect. 3 we present one after the other the proposed model extensions, Sect. 4 includes the complete updated model, and Sect. 5 accounts for a number of numerical simulations exhibiting aspects and properties of the overall updated model.

Previous Hybrid Model
The hybrid model that we present in this paper is an extension of a previously proposed hybrid model detailed in Sfakianakis et al. (2020), Franssen et al. (2021).Here, we provide a brief overview of the model.
The aforementioned hybrid model is comprised of two subsystems: the first incorporates the model components represented by their corresponding physical densities, including ECCs, the etracellular matrix (ECM), and the matrix-degrading enzymes (MMPs).The second subsystem focuses on the mesenchymal-like cancer cells describing them as individual/solitary migrating cells.

Density Subsystem
We consider a Lipschitz set ⊂ R 2 or R 3 and denote the (scalar) densities of the ECCs, ECM, and the matrix-degrading Metalloproteinases (MMPs) by c E (x, t),  v(x, t), m(x, t) where x ∈ respectively.

ECC Density
The ECCs do not migrate actively.It is rather assumed, in this version of the model, that the colony of the ECCs disperses in the surrounding environment due to mechanical forces, e.g.internal pressure due to proliferation.This phenomenon is macroscopically captured through a (small in magnitude) diffusion term.
Furthermore, ECCs can be transformed into Mesenchymal-like Cancer Cells (MCCs) and vice versa by the EMT and MET, respectively.In this version of the model, it is assumed that EMT and MET are random and take place at a fixed rate.It is furthermore assumed that the ECCs proliferate in a logistic fashion, which accounts for the competition for free space among the ECCs, the MCCs and the ECM macro-molecules.These assumptions lead to the formulation of the following equation where c M ≡ c M (x, t) represents the density of the MCCs, is the indicator function of the set S defined as: Here EMT occurs in a randomly chosen subset of the whole domain, E(t) ⊂ , with a fixed rate μ E .In a similar fashion, it is assumed that every (solitary) MCC might undergo MET independently from the others.The MCCs that undergo MET give rise to the subset M(t) ⊂ ; MET is assumed to occur with a fixed rate μ M .Refer to Sect.2.4 for a full discussion of the EMT and MET operators.

MMPs Density
The MMPs are assumed to be produced by both types of cancer cells and to diffuse in the cancer microenvironment.It is also assumed that the MMPs decay at a constant rate.The time evolution of their density is given by the PDE:

ECM Density
The ECM is represented by the density of collagen macro-molecules and it is modelled as a nonuniform and immovable component of the system that neither diffuses nor otherwise translocates.Additionally, it is assumed that the ECM is degraded by the action of the cancer-cell/MMP complex (both epithelial and mesenchymal).Finally, for the sake of model simplicity, it is assumed that the ECM is not reconstructed.The evolutionary equation of the ECM density is the following The assumption that the ECM is degraded by the cancer cell-MMP complex, rather than simply by the action of MMPs alone, is motivated by the fact that cancer cells require MT1-MMP, a specific type of membrane-bound MMPs, for invasion to take place (Farideh et al. 2009).

Solitary Cell Subsystem
The MCCs are represented by a system of solitary cells indexed by p ∈ P = {1, 2, . . ., N }, N ∈ N. Due to EMT and MET, the number N of solitary cells varies in time and, therefore, N = N (t).The MCCs are represented as point masses; x p (t) ∈ R 3 (or R 2 ) represents their position and m p (t) ≥ 0 their masses.Overall the set of solitary cells is described by (5) The migration of MCCs follows a biased random motion comprised of two assumed to be independent processes: a directed migration part that represents the haptotaxis response of MCCs to gradients in the density of the ECM (drift term), and a stochastic component that represents the undirected kinesis of MCCs as they sense their environment, which is understood as Brownian motion (diffusion term).
The Brownian motion assumption is clearly a simplification justified by the random walk-like migration that the cells exhibit, see also Stokes et al. (1991).Based on these assumptions, the migration of the solitary MCCs obeys the following SDE where X p (t) = X p t denotes the position of MCCs where W p t is a Wiener process with independent components.The drift and diffusion coefficients encode the modelling assumptions made about the nature of the directed and random components of the motion.We refer to Sfakianakis et al. (2020) for more details regarding (6).With X p t we denote the stochastic process describing the position of MCCs, satisfying SDEs of the form (6) and with x p we denote the physical position of the solitary cancer cells represented by the set P(t) in (5).

Modelling Reactions of MCCs
It is clear that (6) describes only the migration of the MCCs.Still, these cells participate in a number of reaction processes such as MET, EMT, production of MMPs, degradation of the ECM, and more.
These reaction terms are accounted for as follows: the MCCs undergo MET randomly, after which they are removed from the set P in (5) and transformed to a density, via the solitary-cell-to-density operator (introduced in Sect.2.3), that is then added to the existing density of ECCs.Conversely, parts of the ECC density undergo EMT in a random fashion; this is initially transformed into MCC density, and subsequently, Besides EMT and MET, the full set of solitary MCCs is translated into density, denoted by c M in (1), through the solitary-cell-to-density operator described in Sect.2.3 below, and subsequently employed in the evolution of the other components of the system that are described as densities.This is because the density of MCCs participates in the proliferation of ECCs (1), the production of MMPs (3) and the degradation of the ECM (4), which are all described via the density submodel.

Phase Transitions Between Densities and Solitary Cells
The phase transition between the density description of the macroscopic ECCs and the solitary MCCs is conducted by the solitary-cell-to-density and density-to-solitary-cell operators.To proceed, it is assumed that the domain is regular (e.g. a rectangle in two dimensions) and is sufficiently large to be split into equivalent partition cells M i = i∈I M i . (7)

Density-to-Solitary-Cell Transition
For a given density function c(x, t), the density-to-solitary-cell operator B is given by c(x, t) where at every partition cell M i we assign a cell with mass and position Although the partition cells M i cover the whole domain, they should not be understood as the computational cells of the numerical discretisation.In a typical application of this method, they are set to have the size of the biological cell under investigation.

Solitary-Cell-to-Density
The solitary-cell-to-density operator F is defined as Every solitary cell is assigned to a rectangular domain K p of size K > 0 (typically chosen to be the size of the biological cells), centred at the centre of mass of the cells x p .It is assumed that the mass of the cell is evenly distributed over K p .As K p overlaps with (possibly) several of the partition cells M i , we assign the corresponding portion of the cell mass given by The mean value of c(x, t) over the partition cell M i is denoted by c i (t).The contribution of all solitary cells ( p ∈ P) to the density of partition cell M i as The density function across the full domain is then given by summing c i (t where X M i (x) is the indicator function defined by equation (2).Note that the transition between the two phases is mass conservative in all the various stages.Note furthermore, that there is a high level of versatility built into the two operators and, in effect, in the hybrid description of the cancer system.For example, the operators B and F allow the materialisation of individual cells of various masses or of multiple cells at a specific location, or even of the complete conversion of cell-density profiles into collections of individual cells and vice-versa.In the current paper though, we use these operators primarily for the EMT/MET processes and to materialize individual solitary cells with mass and volume values corresponding to those of the HeLa (Henrietta Lacks) cell line.Refer to Fig. 1 for a graphical representation of these operators in 2D.

The EMT and MET Operators
In this section, we present the mathematical description of the EMT and MET processes which, along with the density-to-solitary cell and solitary-cell-to-density operators, allow for coupling between the two cancer cell phenotypes.

EMT Operator
It is assumed that a randomly chosen portion of the ECC density, c E E MT , undergoes EMT to give rise to a density of MCCs This newly created MCC density (c M E MT ) is then immediately transformed into solitary MCCs via the density-to-solitary-cell operator, described in Sect.2.3 where x M, p , m M, p are the position and mass of the newly created solitary cells and P E MT is the corresponding set of indices.Following this, the collection of existing solitary MCCs is updated with the addition of the newly created solitary MCCs.This is given by the disjoint union where P new is the re-enumeration of the set P P E MT .Overall, the EMT operator reads as

MET Operator
It is assumed that each of the solitary MCCs undergoes MET randomly, to produce a set of solitary ECCs Subsequently, the solitary ECCs are transformed into a density via the solitary-cellto-density operator Overall, this can be written in operator form as (21)

Updates and Extensions to the Previous Model
The model updates that we propose in this paper encompass several key modifications including the incorporation of nonlinear diffusion-porous medium type-in the time evolution equations of the various cell types.Additionally, we address the role of CAF cells and of TGF-β in the tumour microenvironment.Furthermore, we refine the modelling of the EMT and MET operators and of the SDEs describing the migration of the MCCs.Lastly, we broaden the scope of our hybrid model by extending it to a multipleorgan metastasis framework, thereby facilitating a more comprehensive understanding of the progression and dissemination of cancer across the whole organism.

Nonlinear Diffusion
The equation for the ECC evolution in the previous version of the model, (1), employed linear diffusion which (conditionally) leads to infinite propagation speed that the solutions exhibit.This implies that ECCs could spread instantaneously across the tissue, contradicting the actual biological nature of cancer cell migration, especially considering the timescale of cell migration and tumour growth.To address this discrepancy, our revised model introduces a simple nonlinear diffusion term, effectively resolving this issue.
In general terms, the diffusion functions D(u) in the nonlinear diffusion equation du dt = (D(u)∇u) can be categorised into degenerate and non-degenerate.Degenerate diffusion functions are characterised by their degeneracy at zero, i.e., D(0) = 0.The degeneracy ensures that, for any compactly supported initial condition, there is no diffusion at the interface where u = 0.As a result, the compactness of the initial conditions is maintained at all times.On the other hand, non-degenerate diffusion functions, such as linear diffusion, do not possess the property D(0) = 0.As a result, they do not create a sharp interface, and under certain conditions, can lead to infinite propagation speed.For a detailed discussion and comparison between degenerate and non-degenerate diffusion functions, we refer to Vazquez (2007).This model extension, incorporating nonlinear degenerate diffusion functions, was first explored in Williams (2020), where various such functions were examined.Building on the remarks therein, we propose here the use of a Porous Medium type Equation (PME).The PME represents a natural extension of the classical heat/diffusion equation and is given in general form as where n > 1.Assuming no-negative initial conditions, u remains non-negative for all times.Such would be the case in biological settings where u represents cell density, the PME ( 22) can be written as Clearly, the PME is degenerate as D(0) = 0 and admits solutions with a clear interface, compact support, and finite propagation speed.To demonstrate this key difference in behaviour we have plotted in Fig. 2, typical solutions for both the Heat Equation and the PME.The solution in the Heat Equation spreads asymptotically in space, observing, hence, infinite propagation speed.In the contrary, the PME develops a distinct interface, with steep sides that propagate with finite speed.The PME, as its name suggests, was originally used to describe the flows through a porous medium, but has since been used in many other contexts (Vazquez 2007).In particular, porousmedium diffusion has been found to fit well with experimental data of biological cell migration (Sengers et al. 2007;Jin et al. 2016).That is, it removes the problem of infinite propagation speed while introducing minimal additional complexity to the existing model.

ECCs Density with Porous-Medium Diffusion
We introduce the nonlinear degenerate porous-medium type diffusion into the equation for ECC density.To this end, we consider the nonlinear porous-medium type diffusion function to be given by This is the simplest type of porous diffusion function and corresponds to n = 2 in (23).Using a diffusion function of this form, we implicitly assume that diffusion increases linearly with density.This is desirable as we have already stated that we assume ECCs primarily diffuse due to mechanical forces (e.g.pressure) caused by proliferation.
We retain all the original assumptions from (1).However, we replace linear diffusion with the porous-medium type diffusion function given in (24).Additionally, we assume that the ECCs compete with the CAFs which are included in the updated model and are introduced in Sect.3.2.We have also included a haptotaxis term that accounts for the biased motion of ECCs, which is known to be in the direction of increasing ECM density (Sfakianakis et al. 2017;Anderson et al. 2000).Therefore, the updated equation for ECC density is given by The EMT and MET operations are modelled in this current version as random processes through the variables T EMT and T EMT .Moreover, in contrast to the original modelling in Sfakianakis et al. (2020), EMT depends here on the concentration of TGF-β.These updates of the model are discussed in detail in Sect.3.4.

The Role of CAFs
An additional biological component to the invasion-metastasis cascade are the Cancer-Associated Fibroblasts (CAFs).CAFs are known to reconstruct/re-model the ECM, secrete MMPs to the extracellular environment (Weinberg 2014; Xing 2010; Erdogan and Webb 2017), and are one of the main producers of TGF-β, which is responsible for the EMT procedure and will be discussed in extent in Sect.3.3.We assume that the CAFs perform biased random motion, modelled via a nonlinear diffusion and haptotaxis equation, towards areas of lower ECM density.It is furthermore assumed that CAFs proliferate in a logistic manner, that their proliferation is increased by the presence of ECCs, and that they die at a constant rate.Therefore, the evolution equation of the CAFs is given by The following argument justifies the proliferation term: in the absence of ECCs, the CAFs self-proliferate in a logistic volume-filling manner, in which they compete for resources and space with themselves, ECCs, MCCs and ECM.However, we also assume that more CAFs are produced in the presence of ECCs.This accounts for the transdifferentiation of ECCs into CAFs, as well as the recruitment by ECCs of healthy fibroblast cells and their subsequent transformation into CAFs.Therefore, in areas of high ECC density, we would also expect high CAF density.Other models have explicitly included transdifferentiation (Sfakianakis et al. 2017), or healthy (nonactivated) fibroblasts which then become activated to produce CAFs (Smellie 2022).However, for simplicity, it is assumed that these effects are captured by the modification to the proliferation term.In the case of wound healing, fibroblast cells are directed towards areas of low ECM density as they are capable of remodelling the ECM and thus helping to repair the wounded site (Darby et al. 2014).We have assumed that this will remain the case for CAFs and have modelled this behaviour with a haptotaxis term promoting migration in the direction of lower ECM density.

The Role of TGF-Ť
he cellular (re-)programming processes EMT and MET are not random, they are rather initiated through complex interactions involving extracellular signalling molecules.Specifically, the involvement of TGF-β has been recognized as critical in triggering EMT (Katsuno et al. 2013;Kalluri and Weinberg June 2009;Bierie and Moses 2006).In this study, we focus solely on the role of TGF-β in EMT and do not consider the impact of other signalling molecules.While this is a severe simplification of the biological reality, it allows us to investigate the specific role of TGF-β in the context of EMT and MET dynamics.
From a modelling point of view, we make the assumption that TGF-β is primarily produced by CAFs, while neglecting the contribution of ECCs and MCCs.We further posit that TGF-β diffuses freely in the environment, a condition that stands in contrast to the nonlinear diffusion assumption applied to cancer cells.This choice is substantiated by the differing diffusion time scales that molecular components-unlike cellular ones-exhibit within the tumour microenvironment.Additionally, we assume that the decay of TGF-β occurs at a constant rate.
Hence, the governing equation for the spatiotemporal evolution of TGF-β density is given by

The Updated EMT and MET Operators
Another major extension in the model concerns the EMT and MET operators in (25), which play a crucial role in capturing the transitions between epithelial and mesenchymal phenotypes.In the previous iteration of the model, the corresponding operators were significantly simplified.However, weit is understood that EMT and MET are complex processes influenced by a multitude of factors, including cellular heterogeneity, stochasticity, and the presence of signalling molecules.Specifically, in the current paper, we have incorporated the influence of TGF-β, which has been widely implicated in regulating EMT and MET dynamics in cancer.The updated EMT and MET operators account for the stochasticity inherent in these transitions, allowing for a more realistic representation of the phenotypic plasticity observed in cancer cells.

EMT
The biochemical mechanisms governing EMT remain unclear and, accordingly, significant simplifications of the modelling assumptions need to be made.Nevertheless, we assume in the current paper that EMT depends directly on the local concentration of TGF-β and, in addition, that it takes place in a random fashion.Namely, we assume that a minimum concentration of TGF-β is necessary for EMT to occur.Still, this or even higher levels of TGF-β are not sufficient to trigger EMT.To account for this uncertainty, we make the modelling assumption that EMT occurs in a random fashion.Overall, we denote by T EMT the amount of ECCs transitioning to MCCs as a Poisson random variable with parameter λ > 0 where the probability density function of P(λ) is given by In the case of EMT, we consider λ = ζ(b)τ , where ζ : R → R is the rate at which EMT events occur that depends on the concentration b of TGF-β in the environment and τ ∈ R is a small time interval.In other words, T EMT ∼ P(ζ (b)τ ) accounts for the probability that a events will happen in the time interval τ with rate ζ .The rate ζ is described as a shifted logistic function where L, k, b T ≥ 0 are constants and b is the density of TGF-β.We choose the shifted logistic function as a switch mechanism that will ensure the necessary condition for EMT to happen, meaning that whenever b > b T the rate ζ becomes strictly positive and reaches its highest value L. The values of T EMT will fluctuate at different times during the evolution of the system, since the rate ζ takes values in the interval [0, L], giving us either zero when the threshold b T is not exceeded, or a proportion of ECCs that undergo EMT when it does (Fig. 3).We further consider that ECCs tend to undergo EMT only at times when the density of the solid tumour c E is sufficiently large.Hence, we account for this restriction by multiplying T EMT by a characteristic function of the form: The constraint (31), although imposed on the density profile c E , it is translated in a very natural fashion to a constraint over the mass through the underlying partition cells K p and M i that account for the size/volume of the corresponding biological cells, see ( 9) and ( 12).Practically, the partition cells K p and M i either coincide with the computational cell of the corresponding numerical method, or, when coarser, the computational grid is is coarsened to the appropriate size.This coarsening of the computational grid is used to calculate the new density over the coarse grid and the corresponding mass.Note that M i and K p have the size of the biological cells under investigation and, as such, they are coarser than the numerical discretisation of the domain, see Sect.(2.3).

MET
Along with EMT, its inverse cellular programming MET, is considered to be a driver in the dissemination of carcinomas through the formation of metastases in the primary and secondary locations of the organism.The occurrence of MET depends on the type of cancer cells and tissues under investigation, and has become a target of various clinical trials (Wood et al. 2021).Due to the lack, though, of a clear biochemical triggering mechanism for MET, and for the sake of simplicity of presentation, we do not include in the current work any dependence on extracellular cues such as TGFβ.We will, though, assume that every solitary MCC undergoes MET randomly and independently from the others.More specifically we will assume that MET occurs in any given MCC according to the Poisson random process with a fixed rate r > 0 for any time interval τ i.e.
where r > 0 is constant.
For each one of the solitary MCCs in the subsystem (19), a decision is made, by (32), whether it will undergo MET.Then, with the use of the solitary-cell-to-density operator F, the system is updated as shown in (21).We account for this in a similar way as we did for EMT, by multiplying T MET in (25) by a characteristic function as in (31).
where c M is the density formulation of MCCs.

Solitary-Cell Invasion SDEs
We reformulate the solitary-cell submodel so that the SDEs involve more of the biological properties describing the movement of MCCs in the tissue.The ECM plays a crucial role in the migration of solitary MMCs since MCCs move towards regions with a higher density of ECM.This implies that the velocity of MCCs is proportional to the gradient of the density of the ECM.Taking into account the haptotaxis response of the MCCs, and disregarding random effects for the moment, we write the following Ordinary Differential Equation (ODE) for the position, X(t) = X t , of each solitary cancer cell such as: where v denotes the density of the ECM.The drift force now directly depends on the gradient of the ECM through the function μ.The form of this drift force needs to take into account that there is an upper limit in the migration speed of the cells, regardless of the density of the ECM.We denote this maximum cell speed by V thr .This function μ could take the following form: which attains values from [−V thr , V thr ].Furthermore, we reformulate the (deterministic) ODE (34) to an SDE in order to account for the random motion of the cells as follows: where N (t) is a homogeneous Poisson process with rate λ > 0 and {Y i } N (t) i=1 are Independent and Identically Distributed (i.i.d.) random variables.The Poisson process is a counting process that represents the number of independent discrete events that have occurred up to a time t.It is furthermore stationary, i.e. the probability of a random cell turning in a given time interval is the same for all equal-length intervals.Due to these properties, we find that the CPP describes the random turning in the migration of the cells better than the standard Wiener and hence employ it in (36).

Numerical Solution of the SDEs
For the numerical scheme of the new SDE (36), we will use an Euler-Maruyama type scheme that takes the following form: where μ is the drift coefficient given in ( 35) and σ the diffusion coefficient.For every timestep τ > 0 it holds with λ > 0 the rate of the homogeneous Poisson process N (t) defined in (37) and Y 1 a random variable.For more details on the derivation and convergence of the method (38) we refer the reader to Grigoriu (Aug 2009).We note that in every timestep of our numerical scheme, it is necessary to ensure that the cells do not exceed the maximum velocity, denoted as V thr .This situation may arise, for instance, when the combined effects of the drift and diffusion terms in (36) exceed V thr .In cases where V t+τ > V thr, with V t+τ representing the velocity of MCC at time t + τ , we adjust the velocity to comply with the maximum limit.This adjustment is done by rescaling the newly calculated velocity as follows: Here, V new t+τ is the new, rescaled velocity that replaces the original velocity as computed by the scheme.

Multiple Organ Metastatic Framework
In the final proposed extension of the mathematical model, we introduce a complex multiple-organ metastatic framework.The integration of mathematical models in identifying metastases is of significant importance as it can provide a better understanding of the metastatic cascade and potentially offer a predictive approach to complement existing detection technologies.
Within this framework, we conceptualize in this paper an organism as a network of organs, interconnected by a circulatory system that serves as a conduit for cancer cells in their metastatic spread to secondary sites in the organism.
The model encompasses several critical biological phenomena, beginning with the formation of ECC-like tumour in the primary location within an organism.Following the occurrence of EMT, MCC-like cancer cells emerge in the primary organ.These MCCs migrate within the tissue and, with some probability, they infiltrate the circulatory system in a process known as intravasation.We model intravasation by a Poisson process in a similar fashion as in the EMT and MET.
Once in the bloodstream, the cancer cells acquire the designation of Circulating Tumour Cells (CTCs).The hemodynamic environment of the circulatory system imposes significant stresses upon these cells, leading to a substantial rate of destruction among the CTCs.The survival probability of the CTCs is estimated to be 0.1% (Liotta et al. 1976(Liotta et al. , 1977)).In mathematical terms, we model the death rate of the CTCs through a uniform probability distribution.
CTCs that survive the circulatory network, extravasate into one of the downstreamwith respect to the blood flow-organs.The anatomical arrangement and the hemodynamics of the circulatory system play a significant role in determining the sites where cancer cells may extravasate.Still, for the sake of modelling simplicity, we do not account for the blood flow in this paper nor do we limit the extravasation process to organs in the downstream direction of the blood flow.We rather model the extravasation process as a Poisson event, with the choice of the extravasation location following a uniform distribution with respect to the various organs that comprise the virtual organism.
Upon arrest within the secondary location, CTCs typically undergo a phase of dormancy (Phan and Croucher 2020).For the sake of simplicity of presentation though, we do not account for cell dormancy in this paper.We instead opt to treat the newly arrested cancer cells as MCC-like, having in effect the potential to migrate in the tissue or undergo MET.This results in the formation of newly developed ECC-like densities and the engenderment of new tumours; metastasis has occurred.
In addition to the aforementioned simplifications, the geometrical representation of the various organs considered in this paper has been abstractly described as simple 3D polyhedrons of variable density.However, our multiple organ modelling framework could also account for more realistic organ geometries and structures, we refer for instance to Kolbe et al. (2022), Misra et al. (2009), Zrimec and Busayarat (2004).

The Updated Model
In this section, we present the fully updated model by combining all the updates introduced above in Sect.3, and we also incorporate the final two components of the model-the macroscopic equations of MMPs and ECM.These system components have undergone minimal changes in the new model.

Density Formulation
The density equations for the updated ECCs, CAFs and TGF-β have been derived in detail above and are given by ( 25), ( 26) and ( 27) respectively.To complete the updated density submodel, we must also include the density equations for the MMPs and ECM.

MMPs Density
As before, we assume that the MMPs move through the environment via linear diffusion (molecular diffusion) and decay at a constant rate.However, with the addition of CAFs to the model, and because CAFs are known for one of the main producers of MMPs in the TME (Weinberg 2014; Erdogan and Webb 2017), we now assume that CAFs are producing MMPs alongside with ECCs.Thus, the evolution equation for MMPs density is given by with constants D E , ρ F m , λ m ≥ 0. Note that equation ( 41) describes the evolution of soluble MMPs that are secreted in the local extracellular microenvironment and we do not account in (41) for membrane-bound MMPs such as MT1-MMP (Franssen et al. 2021;Sabeh et al. 2009;Egeblad and Werb 2002).Although, the presence of MT1-MMP has been suggested in Sabeh et al. (2009), to be a necessary and sufficient factor for the migration of MCCs.In order to separate the different types of MMPs we will only account for the membrane-bound MT1-MMP, implicitly, through the new degradation term of the ECM, shown in equation ( 42)

ECM Density
For the new density equation of the ECM, we have added a production term with a constant rate that depends on the density of CAFs, since they are responsible for the reconstruction of the ECM (Erdogan and Webb 2017;Sahai et al. 2020).The degradation of the ECM is affected again by the ECC-MMP complex, instead of the MMPs alone, and we have changed the contribution of MCCs due to the effect of the membrane-bound MT1-MMPs.Hence, based on these assumptions, the evolution equation for the ECM density is given by: with constants λ m v , λ e v , ρ F v ≥ 0. Here, K p (t) represents the physical space occupied by the mesenchymal-like cell with index p and X K p (t) is the characteristic function defined in (2).The degradation of the ECM is directly impacted by the MCCs of the solitary cell subsystem rather than the density description of MCCs in (4).

Overall Density Submodel
Combining all the updated equations, the new overall density submodel can be written as a system of PDEs given by

Solitary-Cell Formulation
As in the original model, the MCCs are described via a system of N solitary cells, indexed by p ∈ P = {1, 2, . . ., N }, where N can vary in time.Both the position x p (t) and the mass m p (t) are accounted for and thus the set of solitary cells is again described by (5).However, the evolution of MCCs is now described by the SDE (36), which is introduced in detail in Sect.3.5.

Modelling Reactions of MCCs
Similarly to the original model, the solitary MCCs participate in several reaction processes.These are EMT, MET, proliferation of ECCs and proliferation of CAFs.However, the new solitary-cell submodel still contains no reaction terms.We account for them in the following ways.A certain proportion of ECCs, determined by the local concentration of TGF-β, undergo EMT according to the updated EMT operator introduced in Sect.3.4.This proportion of ECCs is then removed from the ECC density and added to the solitary MCC set via the density-to-solitary-cell operator which remains the same as in Sect.2.3.As stated in Sect.3.4, the MET operator is described by a Poisson process and so a random number of MCCs undergo MET.These solitary cells are then turned into ECC density via the solitary-cell-to-density operator given in Sect.2.3.Furthermore, since the density of MCCs participates in the proliferation term for the ECCs and CAFs, at every time step the entire collection of MCCs must be transformed into a density distribution.

Numerical Experiments
In this section, we conduct numerical simulations to explore the qualitative properties and solution dynamics of the updated model described in Sect. 4. The specifics of the numerical method used for solving the system is not discussed here.For complete details, we refer to Appendix A and B. All computer simulations and visualisations were conducted in MATLAB 2023A [65] with a WINDOWS 11 PRO desktop computer equipped with an 8-core Intel i7-9700K processor clocked at 3.6 Ghz and with 64 GB of RAM.Indicatively, the computational cost for the Experiment 2 discussed later in this section is approximately 22 minutes of computer time per 1 day of biological time.It should be noted that due to the inherent computational parallelisation capabilities of MATLAB no significant additional computational burden was observed when simulating multiple organ conformations such as in Experiment 4. Still, we have observed significant speed-ups, by up to a factor of ×20, when implementing our numerical methods for single organ conformations in compiled programming languages such as C++ (International Organization for Standardization 2024).Before presenting the simulation results, we set appropriate values for all model parameters and establish suitable initial and boundary conditions.We then perform four distinct computational experiments.Through a series of plots, we visually track the system's changes over successive time intervals.These experiments aim to demonstrate how the updated model behaves under various scenarios.However, more comprehensive studies are necessary to fully grasp the impact of the model's new features on the process of cancer invasion.

Parameterization
The parameter values used in our simulations are detailed in Table 1.Whenever possible, we have sourced these values from existing literature to align our simulations with biological data.However, in several instances, it was not feasible to obtain specific values from literature, necessitating our estimation of some parameters.All values are presented in units of days, which aligns with the timescale over which our simulations are conducted.
The computational domain, denoted as , is defined as [−0.05, 0.05] 3 .This represents a cuboid with each side measuring 0.1cm and a total volume of 0.001cm 3 .All components of the density model are confined within this domain.Consequently, we apply homogeneous Neumann boundary conditions, which are specified below: where n is the outward unit normal vector to the boundary of the domain .In the numerical simulations shown in Sect. 5 we have assumed that the diffusion coefficient σ (X p n ) is the same for all cells and remains constant for all times.We have in addition employed a reflective boundary condition, according to which MCCs that escape the domain are returned to their last known position inside the domain.Experiment 1 Haptotaxis Flow: The first experiment we consider focuses solely on the process of EMT and the subsequent haptotaxis of the solitary MCCs.The parameter values which relate to EMT as well as MCC migration can be seen in Table 1.Further to these dynamics, we consider only the diffusion of ECCs, hence the advection and diffusion coefficients of the other model components are set to zero.Moreover, no MET takes place in this experiment, nor do any other reactions or interactions between the model components occur.The remaining initial conditions and experimental setup are described below.
First, we assume that there are no MCCs at the initial time.Second, we assume a spherical tumour of radius 0.01cm, comprised only of ECCs, is located in the centre of the domain.Within the initial tumour, we assign the density of ECCs to be 3 g cm −3 and outside this tumour the density of ECCs is set to zero.This can be written in terms of a characteristic function, see also (2), as where the set S is given by, Third, CAFs are assumed to be present in a randomly selected 30% of the full domain.
Over this subdomain, Y , the density of CAFs is uniformly distributed over the interval [0, 0.001].
Fourth, It is assumed that initially TGF-β is present only in the area surrounding the tumour, i.e., where no ECCs are present.The density is then set according to a uniform distribution over the range [0, 0 Fifth, It is assumed that MMPs are present in the whole domain, and their density is assigned a value based on a uniform distribution over the interval [0, 0.0001] Finally, since the main aim of this experiment is to observe/verify the haptotaxis migration of MCCs, we devise an ECM that is directional with a density gradient increasing towards one of the corners of the domain.Namely, we set Subsequently, the ECM density values are normalised and, lastly, are brought within a biologically realistic range This experimental setup is largely based on an experiment carried out in Sfakianakis et al. ( 2020), which we now extend to three dimensions and introduce the model updates.
From Fig. 4 we observe two key properties of the model.The first is that, as EMT takes place, an initial density of ECCs is transformed into a number of individual solitary MCCs.Secondly, MCCs undergo haptotaxis, which directs their migration into areas of higher ECM density.The updated SDE, has preserved the haptotatic response of MCCs, which is encoded into the model via the advection term in the SDE (36), while the new random component maintains the desired qualitative behaviour of a biased random migration of MCCs.The more detailed modelling of EMT ensures that the process predominantly occurs on the periphery of the tumour.This is because the concentration of TGF-β is initially set to zero within the tumour, and EMT can only occur where the TGF-β threshold is exceeded.

Experiment 2 Cancer Cell Islands
In this experiment, our aim is to reproduce the biologically observed phenomenon of cancer islands forming away from the main body of the primary tumour.Therefore, in this simulation, we consider the full dynamics of the model.All the parameters used in this simulation can be retrieved from Table 1.In particular, unlike the previous experiment, we now consider the process of MET, which occurs with a fixed rate given in Table 1.
The initial conditions for the ECCs, CAFs, TGF-β, and MMPs, remain the same as in Experiment 1 (45)-( 49).However, the initial ECM density is set to randomly vary over the domain according to a procedure developed in Franssen et al. (2021), where we refer the reader for full details.
In the current paper we only provide a brief description of this process and refer to Fig. 5 for a graphical representation.To begin with, an 8 × 8 matrix is created with entries taken from a standard normal distribution, N (0, 1).A number of refinement steps are taken until the resolution of the matrix reaches the desired (computational) resolution of the domain.At each stage, the size of the matrix is doubled to increase the resolution of the ECM.The entries to the new larger matrix are obtained from interpolating the values of the previous smaller matrix with the addition of a small amount of Gaussian noise.Therefore, as the ECM is refined it preserves the initial randomly chosen structure observed in the 8 × 8 matrix, with areas of higher or lower densities appearing in the same regions of the grid no matter what resolution is used.This procedure is extended into three dimensions.However, due to the increased com- putational time that three-dimensional simulations impose, a maximum refinement of 64 × 64 × 64 will be used for the ECM in this experiment.
The simulation results of this experiment are shown in Figs. 6 and 7: The initial concentration of the ECCs is engulfed by the ECM, the CAFs, the TGF-β, and the MMPs.As time progresses, EMT takes place and solitary MCCs are materialised.These MCCs perform a biased random migration, they break free from the main body of the tumour and invade the surrounding tissue.The formation of solitary MCCs via EMT, their overall numbers, and their migration patterns are subject to the overall state of the model and the corresponding parameters.As the phenomenon progresses, the solitary MCCs undergo MET, according to the procedure previously described, acquire epithelial character, and the resulting ECCs start to proliferate.These newly formed ECCs are described in our model through their corresponding cell densities   (Franssen et al. 2021) and give rise to ECC islands outside the main body of the tumour where they continue to grow.

Experiment 3 Growing & Merging Microtumours
In this experiment, we investigate the complex phenomenon of (micro-)tumour merging within a specified environment.The merging process is a multifaceted event that may involve a combination of mechanical interactions, signalling pathways, and modifications in cellular behaviours.It represents a critical phase in tumour development, that potentially leads to more aggressive tumour growth and poses additional challenges for therapeutic intervention (Minchinton and Tannock 2006).
The initial conditions for this experiment are similar to the previous ones.The main difference is that we consider two spheroids tumours rather than one as initial concentrations for the ECCs, namely, where S 1 and S 2 are the following translations of the set S given in (46): No EMT takes place in this experiment, and as no MCCs are present initially, neither does MET.Consequently, the experiment focuses solely on the behavior of ECCs, without the transformation or introduction of new cell states through these biological processes.
The simulation results are shown in Fig. 8, where, as in the previous experiments, we model the primary location of the tumours as a box-like organ.Within this environment, two initially separate tumours are allowed to grow by the processes prescribed by the model ( 43).The growth of these tumours is driven by cellular proliferation.
As these tumors expand, they actively respond to the varying conditions present in the local ECM.This interaction with the ECM's non-uniformities plays a significant role in their development and growth patterns.Over time, as these microtumors continue to grow and adapt to their immediate surroundings, they eventually reach a point where they come into contact with each other.Upon contact, the two initially separate  43).All panels include the ECC concentration (yellow level 0.99 isosurface) seen previously in Fig. 6. a MMPs: as it is their complex with the ECCs that degrades the ECM, we visualise them here in the vicinity of the ECCs.b, c CAFs&TGF-β: visualised over the half x-domain.In both panels, the isosurface reveals that, according to the ICs ( 45)-( 49), the CAFS appear randomly in 1/3 of the whole domain.The same holds for the TGF-β although the dynamics that govern their evolution are more complex (Color figure online) tumors merge.This merging is a critical phase of their development, representing a new stage in the progression of the cancerous growths within the given environment.Following the merging of the two microtumors, they form a single, larger tumor.This unified tumor exhibits grows primarily along its outer surface.This is noteworthy as an indication that a tumour comprised of separate tumour islands grows faster than a single tumour of the same overall size.
This observation implies a key aspect of tumor biology-a tumor that originated from the merging of multiple smaller islands displays a more aggressive growth pattern, particularly at the periphery.This finding indicates that the internal structure and history of a tumor can have a profound impact on its growth trajectory.Understanding these dynamics is crucial for comprehending tumor development and for devising effective strategies for cancer treatment, as it provides insights into how tumors evolve and adapt over time.

Experiment 4 Multiple Organ Metastasis:
In this experiment, we take a further step and expand our mathematical model to account for several interconnected organs within a single virtual organism.This approach allows us to better understand how cancer progresses in the metastatic cascade.
In line with the previous experiments presented in this work, we represent the various organs as cubes of equal size.The ECM is constructed for each organ separately by the process described in Fig. 5 and later employed in Experiments 2 and 3.The various organs are linked though a simplified circulatory network that facilitates the migration of cancer cells between different parts of the organism.While distinguishing between the venous, arterial, and lymphatic systems is crucial in real-world cancer progression, our model simplifies these into a single circulatory network with bidirectional connectivity, implying that cancer cells can move in any direction between interconnected organs.
In their travel through the circulatory system, the cancer cells, now termed Circulatory Tumour Cells (CTCs), are highly likely to be eradicated, with an estimated survival probability of 0.1% in our simulations.
Regarding the initial conditions and boundary conditions, the ECCs and the ECM in the organ where the tumour initially appeared are the same as in Experiment 2, with no initial ECC concentration in the remaining organs.
The simulation results, presented in Fig. 9 clear exhibit the metastatic cascade as is perceived by the this modelling framework.Namely, a tumour grows in a primary location of the organism.EMT takes place leading tot he invasion of MCCs in the local tissue (first row).Gradually the MCCs enter the blood stream and spread to secondary locations in the organism.Both in the primary and secondary sites, the MCCs undergo (potentially) MET giving rise to ECCs (second row).The new formed ECC islands grow due to proliferation and give rise to a number of micrometastases (third row).These micrometastases serve as precursors to larger and and more dangerous metastases which with time merge and form larger tumours (forth row).
As noted in the previous experiment, the simulations presented here reveal again that when cancer cells initially form separate, disjoint islands, their overall volume tends to increase more quickly compared to a scenario where they start as a single, unified mass.This difference in growth rate is not just marginal, indicating a fundamental aspect of tumor biology.Understanding these dynamics is crucial for developing effective therapeutic strategies, as it highlights the importance of considering not just the size but also the distribution and structure of cancerous growths.

Discussion
In this paper, we introduce an updated version of the previously proposed genuinely hybrid multiscale tissue cancer invasion model (Sfakianakis et al. 2020;Franssen et al. 2021).This enhanced model replicates more naturally the transition from the epithelial proliferation strategy original of the ECCs to the individual invasion strategy of the MCCs.The model is comprised of a system of PDEs and SDEs that describe the progression of the ECCs and MCCs, respectively, while incorporating the phenotypic transitions of EMT and MET.
A significant upgrade in our model is the integration of nonlinear degenerate diffusion-similar to a porous medium-in the evolution equation of the ECCs and the other living cell components.This model enhancement contributes to a more realistic representation of the tumour growth in the tissue microenvironment, addressing the issue of infinite propagation speed that is commonly encountered in linear diffusion models.
Another significant extension is the updated modelling of EMT and MET.We have considered the role of TGF-β as an EMT trigger and detailed the modelling process including the inherent randomness associated with MET.Moreover, we have introduced CAF cells into our model to reflect their role in shaping the tumour microenvironment through the reconstruction of the ECM.
We have furthermore, updated the structure of the SDEs that drive the migration of MCCs.Namely, we have accounted for a Compound Poisson Process as the core modelling framework of the stochasticity in the change of direction in the migration of the cells.This approach is in contrast to the typically employed Wiener processes in problems of a similar nature and we have opted for the Compound Poisson Process for its greater biological relevance.
Another substantial extension of the hybrid model introduced in this paper is a multiple-organ metastatic framework.Namely, we have developed a basic circulatory network that links various organs, each being modelled by a "local" version of the hybrid model.This modelling framework allows MCCs to enter the bloodstream and, if they survive the stresses of the circulatory network, be transferred to secondary locations in the organism.Once there, the MCCs may undergo MET, start to proliferate, and give rise to metastatic tumours.
We have validated these model extensions through several numerical experiments presented in Sect. 5.These simulations demonstrate the updated model's capabilities to reproduce phenomena that are biologically relevant and generate predictions consistent with experimental observations.They highlight the model's potential utility in exploring specific facets of cancer biology and treatment response.
At this stage of cancer progression, detecting these micrometastases is crucial in both understanding cancer and in enhancing the efficacy of cancer treatment.Because of their very small size micrometastases (also known as occult metastases) usually cannot be detected by current imaging modalities such as X-ray, CT scan, PET scan, and MRI (Pantel et al. 1999(Pantel et al. , 2009)).However, as noted in Mao (2022), applying AI (artificial intelligence) techniques such as deep learning to data from combined imaging modalities may prove to be much more effective in detecting micrometastases-a so-called AI-based radiomics approach (Afshar et al. 2019;Chen et al. 2021).Such advanced technologies are essential for their identification both preand intra-operatively.The early detection can significantly influence adjuvant therapy decisions and potentially improving patient outcomes, see e.g.Pantel et al. (1999), Smith et al. (2019).Given the difficulty in detecting micrometastases using currently available imaging modalities, researchers are investigating other potential markers which may prove fruitful such as circulating tumour cells (CTC), circulating tumour DNA (ctDNA) and exosomes (Mao 2022).One other emerging area of research which may prove effective is the phenomenon of vascular co-option which has been shown to be important in the survival of metastatic colonies (García-Gómez and Valiente 2020).Research has demonstrated that extravasated metastatic cells in the brain do not immediately migrate away from the capillaries (whence they have emerged) but rather remain in physical contact.This cell-cell connection and interaction has been termed vascular co-option and is believed to take place for micrometastases in other organs as well as the brain.Current research suggests that vascular co-option is very important for the viability of dormant metastatic cells and their transition to macrometastases, so much so that the phenomenon has been proposed as a hallmark for metastasis initiation (Kienast et al. 2010;Er et al. 2018;García-Gómez and Valiente 2020).Vascular co-option is also being considered as a potential source of biomarkers and new targets for micrometastases (Bentolila et al. 2016).It may also be used as a target for therapeutic strategies which aim to block the transition from micrometastases to macrometastases, thereby halting the metastatic spread.Given this information and the fact that the locations of the micrometastases vary greatly over time in the simulation results, a natural extension and development of the work would be to incorporate the vascular structure of the secondary tissue into the model.
Incorporating mathematical models into the detection strategies could provide a predictive framework to identify micrometastases more accurately and more efficiently, thus offering a powerful tool to complement current technologies by guiding surgical and therapeutic decisions towards more personalized cancer care.
We envision further enhancements, particularly the integration of immune system interactions within the tumor microenvironment and the refinement of the circulatory system for enhanced physiological accuracy, hold great promise for advancing theoretical research and potentially clinical applications in cancer treatment.These developments will not just be incremental improvements but significant steps towards transforming our model into a comprehensive tool that can navigate the complexities of cancer biology.By more accurately simulating the interplay between tumors and the body's defenses, as well as the spread of cancer through vascular networks, it is expected that the model will impact clinical decision-making and optimized treatment strategies.These developments represent a promising direction in cancer research, potentially leading to more nuanced insights and improved therapeutic strategies.
Table 2 Butcher tableaux for the explicit (upper) and the implicit (lower) parts of the third order IMEX scheme (B.4), see also Kennedy and Carpenter (2003)  as the numerical scheme we employ is (partially) FV, raising its accuracy to the second order necessitates the use of flux limiters for the interface reconstruction of the numerical fluxes.Out of a large number of flux limiter options, we have found that the Minimized-Central (MC), see Van Leer (1977), constitutes a robust and efficient choice.
Before solving (B.2), we re-organise its terms in implicit and explicit (IMEX splitting) and, accordingly, (B.2) takes the form ∂ t w h = I(w h ) + E(w h ). (B. 3) The actual IMEX splitting depends on the problem at hand, but in a typical case the advection terms A are treated explicitly in time, the diffusion terms D implicitly, and the reaction R terms either explicitly or implicitly, depending on their stiffness.In the problems that we encounter in this paper, all reaction terms have been resolved explicitly in time.
The semi-discrete problem (B.3) is now solved using a diagonally implicit RK method for the implicit part I(w h ), and an explicit RK for the explicit part E(w h ).Altogether, we solve (B.3) using the additive RK scheme where s = 4 are the stages of the IMEX-RK method, E i = E(w i ), I i = I(w i ), i = 1 . . .s, and { b, Â}, {b, A} are the coefficients for the explicit and implicit part of the scheme respectively.These coefficients can be found in the Butcher Tableau in Table 2, cf.Kennedy and Carpenter (2003).As a final stage of this method, we solve the linear system in (B.4) using the Iterative Biconjugate Gradient Stabilised Krylov subspace method, see Krylov (1931), van der Vorst (1992).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 1
Fig. 1 Schematic representation of the relation between density profiles and solitary as used by the operators (8) and (11).Showing on the x − y plane the centre of mass x p of the individual cells along with the rectangular domain K p that represents the space that the cell occupies.The density of the cells is shown as the piecewise constant over K p and their masses m p as the point in the third dimension above x p (Color figure online)

Fig. 2
Fig. 2 Numerical solutions of the a 1D Heat Equation (HE) and b Porous Medium Equations (PME) over time.In c and d are shown the 10 −5 , 10 −10 , 10 −15 isosurfaces (outward order) of the HE and PME in 3D and at the same final time.A relative examination indicates the compactness of the solution's support the PME, justifying thereby the finite propagation observed by it's solution.This becomes more evident by the 3D isosurfaces, with the smallest one, 10 −15 (close to double precision floating-point number accuracy) being distinctively broader than the 10 −10 and 10 −5 ones in the HE but not in the PME solution.Figure source for (a) and (b) (Williams 2020; Harbour 2022) (Color figure online)

Fig. 4
Fig. 4 Simulation results for Experiment 1-Haptotaxis Flow.Shown here is the time evolution of an initial ECC density (yellow isosurface at density value 1) depleted by EMT and the resulting solitary MCCs (red dots) migrating within the ECM (grey background).The ECM has a clear gradient towards the negative x and y and positive z direction (50)-(52).a EMT takes place in an initially spherical ECC tumour giving rise to a number of MCCs that escape from the main body of the tumour.b The newly formed MCCs perform a persistent random motion and migrate within the tissue towards higher ECM densities, i.e. in the direction of decreasing x and y and increasing z. c, d As EMT continues, more and more ECC density is depleted until the initial tumour is completely transformed into solitary MCCs.e Colourbar for the ECM density, common to all subplots (Color figure online)

Fig. 5
Fig. 5 Construction of the (initial condition) ECM density employed in Experiments 2, and 4. The process starts with random values over an 8 × 8 grid; these introduce the basic structure of the ECM.In every stage of the construction process, the grid is bisected with the new values attained by averaging the neighbouring values of the previous coarser grid with the addition of some additive noise.The process stops when the required grid size is reached.At this stage, the values are scaled to the proper biological range.Figure source and reference(Franssen et al. 2021) Fig. 5 Construction of the (initial condition) ECM density employed in Experiments 2, and 4. The process starts with random values over an 8 × 8 grid; these introduce the basic structure of the ECM.In every stage of the construction process, the grid is bisected with the new values attained by averaging the neighbouring values of the previous coarser grid with the addition of some additive noise.The process stops when the required grid size is reached.At this stage, the values are scaled to the proper biological range.Figure source and reference(Franssen et al. 2021)

Fig. 6
Fig. 6 Simulation results for Experiment 2-Cancer Cell Islands.Time evolution of the ECC density (yellow level 0.99 isosurface), the solitary MCCs (red points), and three slices of the ECM (grey in the background).a An initial spherical tumour consists solely of ECCs.b The ECCs undergo EMT to produce MCCs, these invade the local tissue and begin undergoing MET forming the first tumour islands.c, d The number of tumour islands increases as more MCCs undergo MET.The size of these tumour islands also increases, primarily due to the proliferation of ECCs.e Colourbar for the ECM density, shared by all plots (Color figure online)

Fig. 7
Fig. 7 Additional simulation results for Experiment 2-Cancer Cell Islands.Shown here are the final time, t = 0.500 d, concentrations of the continuum components of the system (43).All panels include the ECC concentration (yellow level 0.99 isosurface) seen previously in Fig.6.a MMPs: as it is their complex with the ECCs that degrades the ECM, we visualise them here in the vicinity of the ECCs.b, c CAFs&TGF-β: visualised over the half x-domain.In both panels, the isosurface reveals that, according to the ICs (45)-(49), the CAFS appear randomly in 1/3 of the whole domain.The same holds for the TGF-β although the dynamics that govern their evolution are more complex (Color figure online)

Fig. 8
Fig. 8 Simulation results for Experiment 3-Growing & Merging Microtumours.Time evolution of two ECC densities (yellow level 0.99 isosurface) growing and merging in the absence of EMT and hence MCCs and MET (Color figure online)

Fig. 9
Fig. 9 Simulation results for Experiment 4-Multiple Organ Metastasis.this figure illustrates the progression (across rows) of metastatic growth in four different organs (columns).The various organs are represented cubes whose size and colour schemes are the same as in Fig. 6.At t = 0.002 d: An initial ECC tumour has already developed to a notable size in the first organ (column 1), at which point EMT leads to the formation of MCCs.At t = 0.38 d: These MCCs migrate within the first organ and have beagun spreading, via the circulatory network, to secondary organs (columns 2 to 4).At t = 0.72 d: micro-metastases are observed forming in the secondary organs concurently with the growth of the primary tumour.At t = 1 d: there is a continued expansion of both the initial tumor and the developing micro-metastases, indicating progressive metastatic spread (Color figure online) , j I j + τ n āi,i I i , i = 1 . . .s,

Table 1
Parameter values used to perform simulations of the model