Multiphasic modelling and computation of metastatic lung-cancer cell proliferation and atrophy in brain tissue based on experimental data

Cancer is one of the most serious diseases for human beings, especially when metastases come into play. In the present article, the example of lung-cancer metastases in the brain is used to discuss the basic problem of cancer growth and atrophy as a result of both nutrients and medication. As the brain itself is a soft tissue that is saturated by blood and interstitial fluid, the biomechanical description of the problem is based on the Theory of Porous Media enhanced by the results of medication tests carried out in in-vitro experiments on cancer-cell cultures. Based on theoretical and experimental results, the consideration of proliferation, necrosis and apoptosis of metastatic cancer cells is included in the description by so-called mass-production terms added to the mass balances of the brain skeleton and the interstitial fluid. Furthermore, the mass interaction of nutrients and medical drugs between the solid and the interstitial fluid and its influence on proliferation, necrosis and apoptosis of cancer cells are considered. As a result, the overall model is appropriate for the description of brain tumour treatment combined with stress and deformation induced by cancer growth in the skull.


Introduction
Lung cancer is one of the most frequent cancer types and additionally has the highest death probability for male patients, cf. the report of the Robert-Koch Institute (Kaatsch et al. 2014). One reason for the high mortality results from metastasation, as lung cancer is not only malignant but can also rapidly spread throughout the body. As a consequence, the metastatic spreading to remote organs and the affection of their functions is often tremendous. This motivates the current contribution to model the metastasation process in brain tissue and its related coupled interactions, such as a chemotherapeutic treatment.
Basically, cancer cells emerge from regular healthy cells, where an initial transformation towards a malfunctional cancer cell has occurred caused by specific mutations in the deoxyribonucleic acid (DNA). Usually, malfunctional cells undergo apoptosis, the programmed cell death initiated by specific triggering signals. In contrast to this, the apoptotic behaviour of most cancer cells changes in such a way that there is no reaction to apoptotic signals. Thus, the ability to overcome apoptosis is one of the hallmarks of cancer, cf. Hanahan and Weinberg (2000). In addition to the omitted cell death, cancer cells can acquire properties to successfully spawn towards different parts of the body. During the evolution of a tumour, cell division increases the number of cancer cells until a supply of the growing cell cluster with nutrients and oxygen is no longer possible by diffusive processes alone. Thereafter, the growth of cancer clusters or tumours, respectively, either stagnates or blood vessel growth (angiogenesis) towards the metastases is initiated, cf. Geiger and Peeper (2009). In the latter case, tumours are further supplied via the blood system, thus creating a progressive growth process, while the ability of non-restricted cell division increases the possibility for further mutations in cancer cells. Additionally, cancer cells disseminate from the Dedicated to Gerhard A. Holzapfel on the occasion of his 60th birthday.
tumour and migrate into the surrounding extracellular matrix (ECM) by either cutting through the ECM by mesenchymal motion or by crawling along the fibres (amoeboid motion), cf. Friedl and Wolf (2003). Invasive cancer cells can enter the blood vessel or the lymphatic system and migrate to distant organs forming metastases in their tissue, cf. Shaffrey et al. (2004).
Lung metastases mainly spread to bones, the adrenal gland, the liver and the brain, cf. Nguyen et al. (2009). The brain, as the organ under discussion, exhibits a highly complex structure formed by a network of cells embedded in an ECM with interstitial and cerebrospinal fluid and blood vessels. While the brain has essential functions such as memory and coordination, cf. Trepel (2017), a disturbance of brain tissue by metastases can dramatically reduce the performance of the brain and may result in life-threatening effects. However, a colonisation of the brain requires the passing of the so-called blood-brain barrier (BBB), a specific coating of the endothelial cells to block the migration of macromolecules and cells into the brain tissue, cf. Wilhelm et al. (2013). The BBB consists of very tight junctions between the endothelial cells that cancer cells may overwhelm by either breaking these junctions or by finding gaps between endothelial cells, cf. Chaffer and Weinberg (2011). However, the new tissue environment forces the infiltrated cancer cells to adapt to the changed conditions. As a result, cancer cells can be found in a dormant or in an active state. In this regard, the amount of cancer cells forming a metastasis, and thus being in an active state, is quite low in comparison with the overall amount of circling cancer cells, cf. Kienast et al. (2010).
In case of a metastasation, its initiation is in equivalence with the onset of tumourigenesis, the growth of primary tumours. Cancer cells start to divide, form initial micrometastases, induce angiogenesis and finally grow towards macrometastases, cf. Fig. 1. As a result, the growing metastases displace the surrounding tissue, thus increasing the local pressure and initiating deformation.
In terms of clinical detection processes, the metastases are at this stage large enough to be identified by the use of standard medical imaging devices. Consequently, this time characterises a typical starting point for a medical treatment. In a conventional medical treatment, a surgical removal of the tumour is conducted which is associated by a consecutive intravascular chemotherapy combined with radiotherapy. In cases where this treatment is not feasible, such as deeply in the brain situated tumours, an alternative treatment can be seen in the direct infusion of therapeutic agents via catheters into the brain tissue, cf. Bobo et al. (1994). During this surgery, the extravascular medication enters the brain through the skull, thus bypassing the BBB that typically would have hindered the passage of therapeutics with high molecular weights.
A continuum mechanical description accompanied by a numerical computation of cancer therapies must catch the basic features of the human body, especially of the organ under consideration. As the brain consists of a porous solid matrix composed of a network of cells embedded in the ECM, it must be understood as a typical porous medium that is percolated by interstitial fluid and blood both in separated compartments but interacting in capillary beds between arteries and veins. Following this, a biomechanical description is naturally based on the Theory of Porous Media (TPM), cf. Ehlers (2002Ehlers ( , 2009. When cancer comes into play, the cancer cells are sticking to the porous solid, while nutrients, medical drugs and further components come along with the blood and the interstitial fluid. In the literature, brain tissue is usually described as a very soft viscoelastic solid with heterogeneous microstructure and different properties in tension and compression, cf. Chatelin et al. (2010), Chaudhuri et al. (2015), Goriely et al. (2015), Zhao et al. (2018) and Zhu et al. (2019). As this is not sufficient for the complex structure of the brain, Hakim and Adams (1965), who have been working on the hydrocephalus problem, presented the early hypothesis that the effects occurring in the brain can only be described by the interplay of several brain tissue components. Decades later, their intuition has been followed by other scientists, such as Nagashima et al. (1987) or Smith and Humphrey (2007), who have been working on biphasic brain descriptions based on a porous brain tissue and interstitial fluid. Furthermore, Comellas et al. (2020) explored the porous and viscous responses of brain tissue material in a very recent publication. Apart from brain descriptions, a variety of further biomechanical problems has been based on multicomponent media, such as the general description of tissue growth by Ricken et al. (2007) or Ambrosi et al. (2011), the investigation of avascular and vascular tumour growth by Krause et al. (2012) and Hubbard and Byrne (2013) or the Fig. 1 Metastasation. The hallmarks of the formation of metastases are a the unrestricted proliferation of cancer cells, b the invasion into the blood vessel system and the extravasation at a metastatic site, c the angiogenesis to overcome nutrient limitations and d the reduced apoptosis or necrosis in the cells, cf. Hanahan and Weinberg (2000) mass transport and deformation of multicomponent tumour growth, cf. Shelton (2011), Sciumè et al. (2013) and Faghihi et al. (2020).
The present contribution aims at modelling the main processes of lung cancer metastasation in the brain, particularly the processes of proliferation and apoptosis. Based on earlier work by Wagner (2014), Ehlers and Wagner (2015) and on a very recent article by Ehlers et al. (2021), a model of three components, namely porous brain tissue, interstitial fluid and blood, is set up in the framework of the TPM. This model is enhanced by experimental data that have been taken at the University of Stuttgart by the group of Professor Morrison at the Institute of Cell Biology and Immunology, cf. Stöhr (2018) and Stöhr et al. (2020). This allows for the identification of the most important model parameters for a reliable description of proliferation and apoptosis processes.
Proceeding from the basic model set-up tailored for the description of metastatic processes within brain tissue, the modelling approach is completed by a set of constitutive equations including a data-based parameter selection and its optimisation. As a result, one obtains a strongly coupled system of partial differential equations that has to be solved numerically. In terms of the present procedure, a finite-element-based implementation is derived within the solver Pandas. 1 By this procedure, two-dimensional (2-d) boundary-value problems are computed that correspond to parameter optimisation studies focusing on proliferation and atrophy. Based on these studies, a combined numerical example is discussed exhibiting the process of metastasation from the initiation of cell infiltration over angiogenesis to the treatment. Therein, examples of different drug concentration levels and their effects on the overall cancer-cell mass are shown. Finally, an additional insight into a threedimensional (3-d) problem is presented leading to future research directions. The tensor calculus and notation used in this article are based on Ehlers (2018a).

The TPM in brief
Based on the brain model by Ehlers and Wagner (2015), the continuum mechanical description of growth and atrophy of brain tumour metastases is mainly based on the mass and momentum balances of the model components. These are the healthy or the cancerous brain material, the blood and the interstitial fluid. As the blood primarily exists in arteries and veins, the fluids are assumed to basically exist in different compartments except of the capillary bed, where they interchange oxygen, nutrients and further ingredients. Thus, the present model is composed of Therein, S stands for the porous brain solid, B for the blood and I for the interstitial fluid.
From the above composition, volume fractions n can formally be defined by the fraction of a volume element dv of over the total volume element dv of . Thus, This definition naturally includes the sum of all volume fractions to yield the so-called saturation condition The governing set of balance equations is given by the mass and momentum balances of all basic constituents, solid, blood and interstitial fluid via with the constraining side conditions cf. Ehlers (2002). In the above equations, = n R is the partial density of given as the product of the intrinsic, effective or real density R and the volume fraction n , is the partial Cauchy stress and the gravitation vector. Furthermore, ̂ and ̂ are the so-called density production and direct momentum production terms coupling the mass and momentum balances of solid and pore fluids. Constraints (5) result from the fact that the continuum mechanical system as a whole is a so-called closed system, while the individual components represent open systems, such that they can mutually interact with each other (Ehlers 2002).
Since the TPM provides a continuum mechanical view onto porous media problems, all terms have to be understood as local means representing the local averages of their microscopic counterparts. In particular, the density production represents an increase of partial density driven by chemo-physical processes like phase transformations or chemical reactions, while the direct momentum production exhibits the volumetric average of the local contact forces acting at the local interfaces (pore walls) between the solid and the pore fluids, on the one hand, and between the fluids, (3) ∑ n = 1 . (4) on the other hand. In addition to the above, div ( ⋅ ) is the divergence operator corresponding to the gradient operator grad (⋅) = ( ⋅ )∕ with as the local position vector. Furthermore, ′ and ′′ are the local velocity and acceleration terms of , respectively, while ( ⋅ ) � is the material time derivative following the motion of .
In the present modelling process, all components are assumed to be materially incompressible meaning that their intrinsic densities R remain constant under isothermal conditions assuming a constant temperature of 37 • C , such that energy balances can be ignored. For the solid skeleton, this basically means, although the intrinsic solid density remains constant, that the partial solid density S = n S SR can vary through variations of the solid volume fraction n S and the solid density production ̂S . Generally, materially incompressible constituents with constant R can be described by volume balances instead of mass balances obtained from (4) 1 through division by R . Thus,

The tumour model in detail
Describing tumour-related processes, we make the assumption that the solid skeleton is composed of healthy brain material SB and tumour or metastatic cancer cells ST . It is furthermore assumed that the tumour can either be in an early avascular stage or, after angiogenesis, in a fully supplied stage. In the avascular stage, the tumour is not supplied by blood vessels and receives nutrients only from the surrounding tissue through the interstitial fluid. After angiogenesis, blood vessels have grown towards the tumour and interchange nutrients and further ingredients with the interstitial fluid and directly (6) (n ) � + n div � =n , wheren =̂R .
with the tumour. This complex situation is simplified by the assumption that the interstitial fluid remains the nutrient supplier, while the effect of the blood is taken into consideration through convenient boundary conditions of the respective numerical model. As a result of these assumptions, composition (1) of the model is specified by also cf. Fig. 2. It is seen from (7) that by the use of the above assumption, the blood is only taken as the blood liquid BL , while the interstitial fluid I is described as a mixture of multiple components I with = {L, N, C, V, D} . Thus, I consists of an interstitial fluid solvent IL with nutrients IN , mobile cancer cells IC , vascular endothelial growth factors (VEGF) IV and therapeutic drugs ID as solutes. As a result, the nutrient exchange between the blood and the interstitial fluid is not taken into the model description, such that mass exchanges by mass production terms only occur between the interstitial fluid and the brain solid, such that ̂B ≡ 0.
Given I = ∪ I , the mass balance (4) 1 can be split into mass balances of the individual components of I yielding has been used (Ehlers 2009). Therein, ′ I is the velocity of I , while ( ⋅ ) � I is the material time derivative following I . Fig. 2 Representative elementary volume (REV) with exemplarily displayed microstructure of tumour-affected brain tissue and macroscopic multiphasic and multicomponental modelling approach Furthermore, rebuilding the mass balance (4) 1 for I by summing up (8) 1 over I obviously yields together with the identity ∑ ( I I ) = has been used. In this setting, I defines the diffusion velocity of I with respect to the motion of I .
The partial densities I of the interstitial fluid components I can also be expressed as follows: Therein, n I is the volume fraction of I as n I = ∑ n I is the volume fraction of the interstitial fluid. As a result, s I = n I ∕n I is the saturation of I in I . While I defines the partial density of I with respect to the volume of the overall aggregate, I I with effective density I R characterises the partial density with respect to the volume of I .
As the partial densities I can be summed up to yield I , cf. (8) 2 , one analogously concludes to the effective density IR of the interstitial fluid reading Furthermore, I I can be expressed by its molar concentration c I m given in mol/m 3 and its molar mass M I m given in kg/ mol via cf. Ehlers (2009). As a result of (11) and (12), the mass balance (8) 1 of the species I reduces to the concentration balance through division by the constant molar masses M I m . As the partial densities I I of the solutes I with = {N, C, V, D} included in (12) 1 are generally negligible with respect to the partial density IL I of the interstitial fluid solvents IL , it is justified to use thus also supporting the assumption of constant IR at constant temperature for the materially incompressible interstitial fluid components.
Following the assumption of materially incompressible constituents , the mass balances reduce to volume balances, such that (6)  where = � − � S defines the seepage velocity of , while S = − 0 characterises the solid displacement vector with 0 as the solid location vector in the reference configuration at time t = t 0 , thus including � S = ( S ) � S . Note that the shape of the fluid balances given in (15) 2 is necessary for the numerical treatment of porous media problems, since the solid material is treated by a Lagrangian description on the basis of the solid displacement vector S , while the pore fluids are described in the framework of a modified Eulerian description relatively to the deforming skeleton by the their seepage velocities .
Furthermore, the volume productions n S and n I are coupled via the mass balance side condition (5) 1 , such that where (6) 2 has been used. Finally, the solid volume balance (15) 1 can be integrated analytically to yield compare "Appendix (a)". In the above equation, n S g can be interpreted as the growth-configuration-based solid volume fraction that is related to the current volume fraction n S through the inverse determinant of the solid's deformation gradient S = + Grad S S with as the second-order identity tensor and Grad S = ( ⋅ )∕ 0 . Note in passing that n S g reduces to the standard reference solid volume fraction n S 0 at t 0 , whenever the integral in (17) vanishes, which is trivially fulfilled in case of n S ≡ 0 . Finally, it should be noted that the notion "growth" does not only include proliferation but also atrophy by necrosis and/or apoptosis.
As the solid skeleton is built from brain and metastatic material, it is assumed that These equations underline the fact that at t = t 0 , the solid volume fraction only contains brain material, while metastatic material only comes into play through the existence of n ST after clustering cancer cells have proliferated, such that the number of cells leap over a certain threshold, thus initialising a micrometastasis growing with n ST . Before this 1 3 happens, cancer cells spread into the brain tissue, while their volume fraction in the brain is not measurable.
Nevertheless, it seems necessary to add the following remark to the above volume balances. Although the components of the overall brain material have been taken as being materially incompressible expressed through constant values of R , there is a mass transfer between the brain solid and the interstitial fluid, meaning that cancer cells and further fluid species can be added to the brain skeleton. However, assuming that cancer cells have the same effective density as the basic brain material and that the masses of the species are negligible, the assumption SR = const. is still justified, while the proliferation process is taken into consideration by an intake of brain volume through n S . Concerning the effective densities R of the interstitial fluid components I , it is seen from (7) 3 and (12) 2 that When cancer cells spread into the brain tissue, this process has to be described by concentrations rather than by volume fractions. Following this, the necessary concentration balances of the solutes I can be obtained with the aid of (13) and (15) together with the saturation condition (3): Therein, and in addition to the standard seepage velocities , the relative velocities I are defined through As these velocities formally also have the character of a seepage velocity, they can alternatively be described by the diffusion velocity of I and the seepage velocity of I , such that I = I + I obviously exhibits the diffusion of I in the moving frame of I .
Apart from the complex matter of volume and concentration balances, the model equations are only complete after having added momentum balances. In the present case, it is sufficient to proceed with the momentum balance of the overall model under quasi-static conditions. Following this, the momentum balances (4) 2 of solid, blood and interstitial fluid are added to yield where inertia terms have been neglected. Furthermore, the overall density or the so-called mixture density is defined as (20) with R from (11).
The equation system described above is only appropriate to solve tumour-relevant problems, once it has been closed by convenient constitutive equations for the Cauchy stresses , the direct momentum productions ̂ and the volume productions n or the density productions ̂ , respectively. In contrast to the Cauchy stresses that have to be found for all constituents, direct momentum productions ̂ and the density productions ̂ are constrained through (5) and ̂B ≡ 0 , such that It should furthermore be noted that the solid volume fraction n S can be computed from its initial value n S 0 , the solid deformation gradient S and the constitutive equation for n S , cf. (17). As a result, the porosity n P = n B + n I is obtained from the saturation condition (3) through n P = 1 − n S . Thus, only the porosity and not the blood and interstitial fluid saturations can be obtained from the relations given so far, such that either s B or s I has to be found by an additional constitutive equation.

Constitutive equations
The brain-tumour model described above has to be closed by a set of convenient constitutive equations that includes, on the one hand, the full complexity of the model and is, on the other hand, as simple as possible. These two seemingly contradictory goals can be met coincidentally with the following additional assumptions: During the avascular growth of a brain metastasis, nutrients and further ingredients reach the cell cluster by diffusion through the interstitial fluid. As a result, the interstitial fluid has to be modelled as a mixture of liquid solvent and various solutes governed by (7). In this state of growth, the blood does not play a dominant role as carrier of nutrients. However, after angiogenesis, the blood (24) s = n n P with ∑ s = 1 1 3 takes over as the main supplier of the growing tumour. To avoid an overburdening of the model with two concurring fluid mixtures in the overall pore space, we have refrained from dealing with an additional blood mixture, thus treating the blood through convenient boundary conditions. With this in mind, the model consists of a porous solid S , the blood B and the interstitial fluid mixture I with ∪ I ingredients that split into IL as the solvent and ∪ I as the solutes.

Thermodynamical restrictions
As the constitutive equations have to fulfil the entropy inequality of the overall model, one has to take a view upon the Clausius-Planck inequality basically yielding cf. Eq. (106) of Ehlers (2002) for a constant temperature . Considering the side conditions (23) and the saturation constraint resulting from the time derivative ( ⋅ ) � S of the saturation condition (3) multiplied by a Lagrange multiplier , the entropy inequality for the present tumour model reads Obtaining (27) from (25), not only the saturation constraint has been used but also the fact that the free energy of a (25) mixture is usually given per unit mixture volume through I I and not per unit mixture mass through I . Thus, the mass-specific free energy of the interstitial fluid mixture is expressed as also cf. (11). Furthermore, note that s I = n I ∕n I is usually also addressed as a saturation as s in (24). Here, it typically acts as a volume fraction of the species I in the mixture I with density IR and can thus be computed via Therein, I I is defined through (12) as the partial mixture density of I with effective density I R .
Based on the principle of phase separation, cf. Ehlers (1989), stating that the free energies of immiscible components of the overall aggregate, as on the microscale, do only depend on their own constitutive variables, the solid free energy S of an isotropic elastic solid only depends on the solid deformation through the deformation gradient S . In case that both incompressible fluids would not only be immiscible but also inert without any ingredients, the free energy of the interstitial fluid I would be constant and the free energy B of the blood would be a function of the saturation s B (Ehlers 2009). However, as the interstitial fluid is a fluid mixture, the volume-specific free energy I I is taken into account instead of the mass-specific energy I . Following this, I I = Inserting these results in (27) leads to the final version of the entropy inequality: Based on (33), the exploitation of the entropy principle yields equilibrium and non-equilibrium solutions that can be obtained by the use of the standard Coleman-Noll procedure (Coleman and Noll 1963), also compare Hassanizadeh (1986), Schreyer Bennethum et al. (2000) and Araujo and McElwain (2005). Before carrying out these results, let us have a look at the mixture species I and the classical relations between free energies I I , molar (mole-specific) chemical potentials I m and osmotic pressures I (Ehlers 2009): Thermodynamic equilibrium: Considering (33) in thermodynamic equilibrium, the terms in brackets are assumed to be independent of their multipliers, the velocity gradients and seepage velocities. Furthermore, mass productions are considered, at the first glance, as non-equilibrium terms.
Given the above, one concludes from line 7 of (33), where in case of thermodynamic equilibrium the term in brackets has to vanish for arbitrary I , that and, as a result, together with (34) have been used, also cf. (29). As the pressure governing I has to equal the partial liquid pressure p I = n I p IR , one easily concludes to As p IR is the effective interstitial fluid pressure, it is easily concluded from (37) that p IR splits into a mechanical part, p IR m , thus also defining the Lagrangian multiplier , and a chemical contribution, the osmotic pressure I .
With this in mind, the equilibrium solution for B yields such that with p dif as the difference pressure between the blood and the interstitial fluid. In addition to the fluid stresses, the solid stress reads where (35)-(38) have been used. Note in passing that p FR is the effective overall fluid pressure, also known as pore pressure, and that the above relation defining p FR recovers Dalton's law (Dalton 1802). It should furthermore be noted that S + n S p FR is the so-called effective stress S ef f defined as that part of the solid stress that governs the stress-strain relation of the solid (Ehlers 2018b), such that the solid stress S reads Finally, summing up the stresses of solid, blood and interstitial fluid given above, one obtains the total stress of the overall model as In addition to the stresses, the direct momentum productions yield in thermodynamic equilibrium where (34)-(38) have been used. Dissipation mechanism: While considering (33) in thermodynamic equilibrium, the terms in brackets have to vanish for arbitrary velocity gradients and seepage velocities. However, beyond thermodynamic equilibrium, the dissipation mechanism allows the terms in brackets to depend on Taking a closer look at the dissipative parts of the momentum productions included in D̂p , one easily concludes to where BS , ILS , IL and I S are positive definite tensors describing a sort of friction between the individual components. Given (46), the dissipation mechanism D̂p yields While the friction I S between the interstitial fluid species and the solid is negligible and has thus been neglected, the remainder of friction tensors reads Therein, BR and ILR are the effective dynamic viscosities of the blood and the interstitial fluid solvent, SB and SIL are the intrinsic solid permeability tensors measured in m 2 with respect to the blood and interstitial fluid compartments. Finally, I is the diffusion tensor of the solutes I in the solution I , while R is the universal gas constant. Note that in case of homogeneous permeabilities and diffusivities, SB , SIL and I reduce to where the intrinsic permeabilities can be substituted by hydraulic conductivities through the relations In the above relation, K S is given in m 2 and k in m/s. Concerning the second dissipation mechanism D̂ , one has to guarantee that D̂ is always positive or zero yielding where the mass-specific kinetic seepage energies of the liquid solvents have been neglected. Furthermore, ̂I L dropped out as the liquid solvent does not contribute to the mass production process. By the use of the mass-specific chemical potentials defines the partial pressure of I as a component of I , the dissipation mechanism D̂ can finally be obtained as Note that the chemical potentials ̄I are related to their molar counterparts in (34) .

Flow and transport of pore-liquid components
Interstitial fluid mixture: As the interstitial fluid I consists of the liquid solvent IL and solutes I , the chemical potentials, osmotic pressures and mechanical pressures of these components have to be specified. Given (34), chemical potentials and osmotic pressures of the species I depend on the volume-specific free energies I I that can be given as with I 0 m as the so-called reference or standard-state chemical potential, such that Based on (37), the effective interstitial fluid pressure is In (37) as well as in (56), there has no difference been made between the liquid solvent and the solutes concerning the definition of the potential I I . As a result, p IR has to be found from the boundary-value problem, while = p IR m results from (56).
Furthermore, combining (46) and (48), the direct momentum productions read Given the above equations, (57) 2,3 result in where IL ≈ I has been used.
Inserting ̂ I together with I from (37) in the fluid momentum balance (4) 2 yields under quasi-static conditions with n I I as the filter velocity of the interstitial fluid. For the solutes I , the same procedure yields with the aid of ̂ I and I from (35) (57) where the diffusion velocity I = I − I together with I ≈ IL has been used. Solving (60) with respect to I leads to In this equation, the extension can be neglected with respect to the fact that I S has been neglected in (46) 3 , cf. the remark between (47) and (48). As I S would basically have the same shape as the friction tensor ILS from (48) 2 , neglecting the extension means to neglect the seepage velocity I of the solutes I compared to the other terms in (61). Formally, I would be governed by Following this, the momentum balance of the solutes reduces to where (55) 2 has been used.
Blood plasma: The momentum balance of the blood has to be combined with the momentum production ̂ B from (46) 1 and the stress tensor B from (38). Thus, While the pressures p BR and p IR have to be found from the boundary-value problem, the blood saturation s B has to be constructed constitutively under consideration of (39). Alternatively, one pressure and s B can be obtained from the boundary-value problem, such that (39) can be used directly for the determination of the missing pressure.
Blood saturation and angiogenesis: As has been said before, the volume fraction n B = s B n P of the blood plasma cannot be found from the solid deformation, only the porosity n P = s I + s B is obtained as a function of det S through n P = 1 − n S with n S after (17). With this in mind, s B can be coupled to the angiogenesis in case that relation (39) between the difference pressure p dif and s B is satisfied and the potential BR B is positive in the range 0 ≤ s B ≤ 1.
Following this, the saturation function should account for blood vessel growth by angiogenesis induced by vascular endothelial growth factors (VEGF) during the cancer cell proliferation as well as for elastic deformations of arterial walls during a change in the local pressure conditions. (61) In this regard, the saturation curve s B depends on the pressure difference p dif , cf. Ehlers and Wagner (2015), and a mass-specific Helmholtz free energy B, angio related to the angiogenesis process. Basically, growth processes can be described by smooth Heaviside or sigmoid functions, cf. Hubbard and Byrne (2013), here defined as For arbitrary but growing x, sig (x) has values between zero and one with sig (x = 0) = 0.5 . By the introduction of the constants , and , the use of the argument x instead of x changes the steepness of the curve with increasing gradients for > 1 and decreasing gradients for 0 < < 1 . Substituting 1 by > 1 lowers the value of sig (x = 0) with increasing values for 0 < < 1 . Finally, substituting exp (x) by exp ( x + ) shifts the curve along the x-axis. Extending (65) in this sense yields Applying the structure of (66) to s B (p dif ) results in where has been substituted by B , by exp ( B B ) and by A B, angio . As p dif is given in Pa, B has the unit 1∕Pa . B, angio is the stored angiogenesis energy per unit mass, such that its unit is J/kg. Thus, A is given in kg∕J . Finally, as B is given in 1/Pa, B has the unit Pa of a pressure.
Initially, the volume fractions of the blood plasma and the interstitial fluid are chosen as n B 0 = 0.05 and n I 0 = 0.20 , cf. Nicholson (2001) and citations therein, such that the initial blood saturation results in s B 0 = n B 0 ∕n P 0 = 0.20 . This value is also obtained by the use of (67) with B = 1863.08 Pa , A = 0.0035 kg/J and B = 0.0035 1∕Pa , thus also establishing an initial pressure difference of The above pressure difference p dif 0 originates from the blood and interstitial fluid pressures in the initial state, where the latter corresponds to p IR 0 = 258 Pa ≈ 1.9 mm Hg and originates from the mean interstitial fluid pressure in the brain, cf. Boucher et al. (1997). Concerning the blood pressure, a homogenisation approach can be taken into account favouring the blood pressure of small blood vessels that are equally distributed within the brain tissue, where larger arteries and veins are neglected. As a result, the implemented blood pressure of p BR 0 = 1725 Pa ≈ 13 mm Hg is below the observed human mean arterial .
blood pressure but accepted as the mean blood pressure in the brain.
Given (67), this equation can be solved with respect to p dif (s B ) yielding the inverse sigmoid function displayed in Fig. 3b. The potential corresponding to p dif is obtained from (39) through integration as In the above potential function BR B , the term in the last brackets cannot formally be integrated, as it is a so-called integral logarithm of the dilogarithm type Li 2 (s B ) . Instead, Li 2 (s B ) can be computed as a sum via The potential BR B is plotted in Fig. 3c.
When the growing cancer cell cluster has consumed most of the nutrients that could reach the cell cluster by diffusion, the cancer cells start to send out VEGF proteins to initialise blood vessel growth necessary for their further proliferation, cf. Finley and Popel (2013). In particular, blood vessel growth triggers the creation of new vessels from which nutrients can diffuse into the interstitial fluid or directly into the cancerous tissue. In this regard, the angiogenesis Helmholtz free energy term B, angio is chosen to be proportional to the VEGF concentration c IV m and the threshold c IV m , thus initiating the blood growth via Therein, the parameter angio = 8.625 × 10 6 (mol J)∕(kg) 2 relates the molar concentration of VEGF to B, angio , while c IV m = 2.5 × 10 −11 mol∕m 3 , both triggering the increase or decrease in the blood saturation s B .

Brain skeleton and growth-dependent solid elasticity
Although brain material generally exhibits a viscoelastic material response, we refrain from including the viscous effects of the brain skeleton with respect to the extremely slow deformations resulting from tumour growth. On the other hand, the viscosity of the interstitial fluid and the blood is fully considered. Thus, the material law for the brain skeleton is described by an elasticity law of neo-Hookean type under consideration of growth phenomena and has to be formulated for the solid stress tensor obtained from (41): Therein, the deformation gradient S can be split, as in thermoelasticity or in elasto-plasticity, into two parts, a purely mechanical part Sm and a growth-dependent part Sg , also compare "Appendix (b)". This split can be motivated on  (17)  reading where n S g =∶ n S 0 det Sg is taken as an a-priori constitutive equation for det Sg . As the integration of the solid volume balance (15) it is evident that det S rules the relation between the solid volume fraction n S of the current configuration and the growth-depending volume fraction n S g . This implies that n S g as a function of growth apparently substitutes the reference volume fraction n S 0 . Only in case that no growth appears, such that n S g = n S 0 with det Sg = 1 , the standard relation det S = n S 0 ∕n S is recovered for materially incompressible solid skeletons. As a result of the above, the multiplicative split of growth-depending materials results in where the Jacobian determinants can be expressed by solid volume fractions via also cf. Fig. 4. As can be seen from this figure, the purely mechanical deformation included in Sm takes place between the reference configuration S 0 (t 0 ) and the current configuration S (t) , while the total deformation represented by S is between the growth configuration S g (t) and S (t) . This is different from the deformation split in thermoelasticity and in elasto-plasticity, where an intermediate configuration S int (t) is situated between S 0 and S and where the purely mechanical, respectively, the purely elastic deformation takes place between S int and S . This is also different from the assumption made by Rodrigues et al. (1994), Ambrosi and Preziosi (2002) or Garikipati et al. (2004) who all assumed the standard decomposition of S in mechanical and growth-dependent parts.
On the other hand, Fig. 4 can also be interpreted differently. In this case, S g can be taken as a time-depending reference configuration, while S 0 would act as an intermediate configuration at t 0 . Anyway, the order of S shifting from S g to S 0 through Sg and from S 0 to S through Sm would then be fully classical.
Describing tumour growth, we make the assumption that Sg is a purely volumetric deformation yielding such that With (75), we proceed from the same type of constitutive equation as has been used in thermoelasticity by Lu and Pister (1975).
Based on (73), (75) and (76), the basic deformation tensors read Note in passing that S and S are the right and left Cauchy deformation tensors situated in g and , while Sm and Sg both correspond to in 0 .
As in thermoelasticity, the strain energy function does only depend on the elastic contribution of the deformation. In thermoelasticity, this is based on the fact that an unconfined body extends volumetrically when heated without exhibiting stresses. However, under fully confined boundary conditions meaning that no volumetric extension is possible, the overall volumetric strain is zero but splits into a thermal and a mechanical part of equal size. Here, the mechanical part induces stresses to compensate the heat extension.
Comparing volumetric growth to extension under heat makes apparent that growth only induces stresses when a Sm .
(77) free growing is hindered. Following this, the solid strain energy function is only based on the elastic part of the deformation governed by Sm . A detailed derivation of the following elasticity law is found in "Appendix (c)". Based on (71), the effective stress can be modified towards with W S ∶= S 0 S as the hyperelastic solid strain energy including the partial solid density S 0 = n S 0 SR of the solid's reference configuration at t 0 . Following the argumentation by Ehlers and Eipper (1999), a nonlinear strain energy function of neo-Hookean type for porous solid materials formulated in the first and third invariant of Sm , I m = Sm ⋅ and III m = det Sm , is given by 0 and S 0 are the Lamé constants defined with respect to the solid reference configuration S 0 . Based on (78) and (79), the solid effective stress reads with Sm = 1 2 ( Sm − ) as the mechanical part of the Karni-Reiner strain S = 1 2 ( S − ) . cf. Ehlers (2018b). In the literature, the neo-Hookean term 2 S 0 Sm is often modified, for example, for the inclusion of anisotropy, cf. Ehlers and Wagner (2015), or for the inclusion of different material properties under tension and compression, cf. Comellas et al. (2020). In the present case, these features are of minor importance, as tumour growth and atrophy initiate pure volumetric deformations.
As the Karni-Reiner strain Sm and the Jacobian J Sm only represent the mechanical part of deformation and strain, these terms have to be substituted by the total deformation represented by S and J S = det S and the growthdependent deformation through (73) and (75). As a result, the growth-dependent elasticity law yields

Basic setting
The mass production terms of the model govern the biological processes occurring during the evolution of tumours. These quantities have to be derived constitutively in a thermodynamically consistent way meaning that they have to fulfil the thermodynamical restrictions (53).
For the further understanding of the mass interaction process, recall that ̂S +̂I = 0 and that ̂I splits into contributions ̂I of the interstitial fluid solutes I with = {N, C, V, D} , nutrients, cancer cells, VEGF and drugs. As the individual ̂I can interchange with the solid skeleton, each ̂I has a solid counterpart ̂S in ̂S , such that Thus, (53) can be rewritten to yield Furthermore, the individual production terms generally split in their gains ̂I ⊕ and losses ̂I ⊖ following Krause et al. (2012): Therein, * is a placeholder for the source of gains and losses. In the following, mass interactions for biologically induced metastatic processes, proliferation, necrosis and apoptosis, cf. Fig. 5, are stated for the solid skeleton and the interstitial fluid.

Mass production terms of the metastatic solid skeleton
Basically, nutrients supply the healthy brain as well as the metastatic tumours with energy for basal reactions as well as for their proliferation. Here, the metastatic proliferation ̂S T ⊕, IN as a part of ̂S corresponds to its nutrient consumption ̂I N ⊖, ST as a part of ̂I.
In the present study, only those nutrients are taken into consideration that are consumed by the cancer cells during proliferation or basal reactions, cf. Guppy et al. (2002). All other substances needed for basal reactions and cell proliferation are assumed to be available in a sufficient amount. In malnutrition states, meaning an undersupply of nutrients, the solid skeleton (cancer and regular cells) will undergo a necrotic process via ̂S T ⊖, IN and ̂S B ⊖, IN . Right before this incident, the cancer cells start synthesising VEGF to initialise angiogenesis. In addition to this, the apoptosis reaction ̂S T

⊖, ID
based on the effect of the therapeutic agent is also considered resulting in a mass-loss term. Following this, the mass production term of the solid skeleton is found as where only ̂S B∕ST IN exhibits gains and losses. Furthermore, all gain and loss processes are interacting with corresponding gains and losses of the interstitial fluid.

Mass production terms of the interstitial fluid
Nutrients: The change in nutrient mass is mainly governed by its consumption from metastatic and freely floating cancer cells through ̂I N ⊖, ST and ̂I N ⊖, IC as well as from healthy brain tissue through ̂I N ⊖, SB . The only increase is due to the triggering angiogenesis based on the VEGF concentration ̂I N ⊕, IV . Therefore, the overall nutrient mass changes via Cancer cells: Like the metastases, freely moving cancer cells proliferate through ̂I C ⊕, IN due to the consumption of nutrients via ̂I N ⊖, IC . On the other hand, the lack of nutrients reduces the cancer cell mass by ̂I C ⊖, IN . Additionally, the interplay with the therapeutic agent leads to a mass loss governed by ̂I C ⊖, ID . Summing up, the mass production of the cancer cells results in VEGF: The increase in nutrient consumption decreases the amount of nutrients in the interstitial fluid. This leads, firstly, to a stop in cancer cell proliferation and, secondly, to necrosis. Counteracting necrosis, VEGF is synthesised within the metastases and released into the interstitial fluid resulting in ̂I V ⊕, ST . Here, VEGF binds to the respective ligand side at the endothelial cells, thus triggering its growth in the direction of decreasing VEGF, cf. Weinberg (2011), Geiger andPeeper (2009) and Potente et al. (2011). However, the synthesis is much higher as its amount is consumed by the growth of endothelial cells with the result that there is a measurable amount of VEGF increase in the interstitial fluid and the blood circulation, cf. Werther et al. (2002) and Kut et al. (2007). This leads to Therapeutic agent: Metastases can only be discovered if their mass, respectively, their size, is above a certain threshold. Upon detection, a drug treatment can be initialised triggering apoptosis of cancer cells in the interstitial fluid through ̂I D ⊖, IC and of metastases in the solid skeleton via ̂I D ⊖, ST . Furthermore, the therapeutic agent is reduced over time through ̂I D ⊖, IL resulting in a half life of the drug of roughly 5 hours in case of a TRAIL 2 derivate within mice, cf. Walczak et al. (1999). In the present study, the half life is considered to 5 h to compensate for the use and clearance of drugs as well as for its depletion into the blood vessels. Thus, The process of angiogenesis couples the mass exchange between the interstitial fluid and the blood which, in contrast to density productions, is described within the constitutive formulation of the saturation function

Results of the metastatic mechanism
Based on the above results, the individual production terms can be grouped together with their respective counterparts yielding where (84) (90) as ̂I N IV is based on the angiogenesis process and ̂I D IL is constantly depleting. Furthermore, it has been assumed that in addition to the general constraint of the mass productions, ̂S +̂I = 0 , the particular mass productions together with their counterparts vanish. Based on this assumption, (83) and (90)  are only related to their specific chemical potentials ̄I N and ̄I D , respectively, such that From a biological point of view, the nutrients would be released from the blood into the tissue of brain solid and interstitial fluid, cf. Abi-Saab et al. (2002). However, the model does not represent the blood as a mixture including interchangeable nutrients that would result in a dissipation relation similar to (91). Thus, the solid chemical potential is omitted in the dissipation inequality for the nutrient gain as it is not the primary interacting term. As a result, ̄I N is assumed negative, thus allowing for a simplification of the model that ensures the biological process and allows, at the same time, for a positive constitutive ansatz function for ̂I N IV . Similarly, the drug loss ̂I D IL is related to metabolisation as well as to depletion into the blood vessels. In the tissue, the metabolisation always results in smaller molecules increasing the dissipation. As a result, ̄I D is assumed negative to allow for a positive ansatz function to equalise the negative sign of the mass reduction. Therewith, the dissipation inequality (83) is not only satisfied in a sufficient manner, but, at the same time, one gains a switch for the biological processes behind the particular productions. No matter, if the productions included in (91) are gains or losses, the corresponding process is only active when the particular inequality is fulfilled. For the gains and losses of the healthy brain material SB , for example, one obtains in dependence of nutrients IN For the gains and losses of all further mass productions, equivalent inequalities are obtained. This implies a threshold for the process and furthermore allows for a positive constitutive ansatz for the individual mass productions. Thereby, it is unnecessary to distinguish between gains and losses, as in case of losses, the order of chemical potentials in (91) changes.

Constitutive setting of mass production terms
Solid skeleton: For the description of the individual mass production terms of the solid, Monod-like kinetics are used for the proliferation in accordance with Ambrosi and Preziosi (2002), Shelton (2011) and Sciumè et al. (2013), whereas atrophy follows a linear behaviour. Following this, the growth and atrophy terms describing reactions to nutrient supply and drug medication can be modelled via Therein, ⊕ ST, max , ⊖ ST, max , ID, max and ⊖ SB, max are maximum proliferation rates with the dimension 1∕s . ̃I D, max is given in m 3 ∕(mol ⋅ s) and characterises the maximum reaction rate. Furthermore, K gr indicates the concentration of c IN m , at which half of the maximum proliferation rate ⊕ ST, max is reached. Moreover, the mass production terms are only taken into account if their considered process is relevant. In case of the proliferation process, this appears when the cancer cells need to have a sufficient nutrient concentration As a result, one obtains for the growth and atrophy of cancer cells in analogy to (94) Kallinowski et al. (1988), Vaupel et al. (1989), Choi et al. (2002) or Cherk et al. (2006). Moreover, N IC is the number of cells per unit mass, given in cells∕kg . In (96) 3 , f proli is a constant factor multiplying the sum of production terms. The remaining term within the nutrient production ̂I N corresponds to the angiogenesis and the nutrient increase of the interstitial fluid. This term reads Therein, the nutrient increase is proportional to the maximum nutrient release rate ⊕ IV, max corresponding to the blood volume fraction. Together with the nutrient concentration c IN m , its initial value c IN 0 m and the molar mass M IN m , it is responsible for the angiogenesis resulting from n B − n B 0 . Recapitulating, the increase in VEGF has been related to an increase in the blood saturation or the blood volume fraction, respectively, and finally to the nutrient concentration. Thus, the metastatic brain tumour triggers the angiogenesis process by the production and release of VEGF into the interstitial fluid. To simplify this process on the modelling side, it is assumed that the cancer cells generate VEGF according to the following equation: The VEGF synthesis is initiated for nutrient concentrations below c IN m and is proportional to the molar mass M IV m of VEGF and the maximum rate ⊕ IV, max given in mol∕(m 3 ⋅ s).
However, the apoptosis is only activated if the threshold concentration c ID m of therapeutic drugs is reached and the Heaviside function H 4 (c ID m ) switches from 0 to 1. This enables a tumour and cancer cell decline due to drugs, while the apoptosis reduces the amount of the therapeutic agents at the same time. As the drug binds to the receptor side on the cell surface, it does not contribute to any further apoptosis reactions. Following this, the mass of the therapeutic agent is reduced as follows: Therein, the two terms of (89) given in (94) 3 and (95) 3 have been summarised to yield ̂I D ⊖, ST∕IC . Furthermore, the dimensionless factor ̆I D, max relates the mass loss of the cancer cells to the therapeutic agent, and ̄I D corresponds to the half life of the drug.

The micrometastatic switch
During the cancer cell proliferation and the infiltration of cancer cells into the brain tissue, a volumetric growth occurs over time and can be expressed as a volume fraction after the growth has exceeded a certain critical cancer cell concentration. This situation characterises the so-called micrometastatic switch. As the cancer cell concentration is inherently not measurable by volume, there exists a related volume given by the volume fraction n IC mainly as a function of the cancer cell concentration c IC m . In order to relate concentrations to volume fractions, recall that partial densities can be expressed by volume fractions and effective densities, on the one hand, and by molar concentrations and molar masses, on the other hand. For the species I of the interstitial fluid, this yields where (10) together with (12) 1 have been used.
In the avascular stage, cancer cells reach the brain as components of the interstitial fluid. Applying this to the general relation (100) 2 yields the volume fraction of metastatic cancer cells in the interstitial fluid of the brain as During proliferation, the number of cells increases until a micrometastasis evolves, cf. Geiger and Peeper (2009). This situation is known as the micrometastatic switch defined by the critical cancer cell concentration c IC m . At this stage, the micrometastasis that has been floating in the interstitial fluid adheres to the solid skeleton, such that it is, on the one hand, . a loss for the interstitial fluid and, on the other hand, a gain for the solid skeleton. Thus, with ICR = STR = SR . Typically, micrometastases are characterised by a specific diameter of approximately 2 mm which can be related to a number of roughly 1.3 million cells, cf. Hermanek et al. (1999). This yields a cancer cell concentration threshold of c IC m = 1.7 × 10 −16 mol/m 3 . Thus, once this concentration is exceeded through a micrometastasis has emerged.
Following this, the micrometastatic switch from the threshold cancer cell concentration c IC m to the initiation of the metastatic volume fraction n ST 0 at t = t ST ≥ t 0 is governed by Further growth is contributed to the metastatic volume fraction n ST via n ST .

Weak form of the governing equations describing the finite-element analysis of growth and atrophy of brain metastases
Targeting the use of the finite-element analysis (FEA) for the numerical computation of proliferation and atrophy processes, the basic governing equations have to be transformed into their weak forms. As the problem is generally highly coupled through the simultaneous action of the overall momentum balance, the fluid mass balances and the concentration balances of the interstitial fluid species, a fully coupled algorithm is strove for, where the system of weak equations is solved simultaneously within the Bubnov-Galerkin method. Starting with the overall momentum balance (21) tested by S , the corresponding weak form reads with the aid of (23) where ̄ = ∑ is the total external load vector acting at the Neumann boundary with outward-oriented unit normal vector , while and Γ are the domain and the surface of the solid skeleton under study.
On the basis of (15) 2 , the volume balances of the blood and the overall interstitial fluid are tested by p R and yield in their weak forms where v = n ⋅ are the volumetric effluxes of the pore fluids .
Finally, the relevant dissolved species I of the interstitial fluid mixture, cancer cells, nutrients and drugs, are described by the weak formulations of their concentration balances (20) tested by c I m yielding Therein, ̄I = n I c I m I ⋅ is the molar efflux of I across the Neumann boundary. As the problem under study is growth-dependent, such that n S cannot be found from n S 0 and det S alone, an additional equation is necessary to obtain n S . This term can either be found by combination of (6) and (17) yielding or by the use of the weak form of the solid volume balance (6) tested by n S , such that Therein, v S = n S ( S ) � S ⋅ = 0 . This means that in a Lagrangian setting, no solid volume flux across the surface of S is possible.
Commenting on the choice of using either (108) or (109) in a numerical study, one must be aware that (108) can only be taken when ̂S is given in an integrable form. Otherwise, one has to proceed with (109) instead.
The above equations are sufficient to compute the primary variables of the model. These are the solid displacement S , the effective pressures p BR and p IR of the blood and the interstitial fluid mixture and the relevant concentrations c I m of the (106) considered interstitial fluid species, while the solid volume fraction, as a result of the integrability of ̂S , is based on (108) and is therefore a secondary variable. The remaining secondary variables, such as the stresses and the momentum productions, can be obtained as functions of these terms, cf. Ehlers (2009) and Ehlers and Wagner (2015). The spatial solution of the overall model makes use of the solver PANDAS and is based on the extended Taylor-Hood elements with quadratic approximation functions for the solid displacement S and linear approximation functions for the remainder of primary variables. In the time domain, we proceed from an implicit Euler scheme.

Parameter identification and optimisation
Due to the large number of simulation input parameters in the governing balance equations of the overall model, a great uncertainty is found for the quantitative description of the simulated processes. This problem can partially be overcome with the aid of raw data directly obtained from experimentalists, and a subsequent identification and optimisation of the most relevant model parameters. In this procedure, a first step is to show what parameters are sensitive for the currently considered process, while, in a second step, the most crucial parameters are optimised. In this study, we rely on cancer cell proliferation and on apoptosis-triggering medication experiments that have been conducted in vitro on cancer cell cultures.
The following parameter identification process is based on the parameter sampling method originating from Morris (1991) together with sensitivity measures introduced by Campolongo et al. (2007). In this context, the deviation between the model and the data is denoted as the model output y. The dependency of y on each parameter is evaluated using the elementary effects method introduced by Morris (1991) and improved in Campolongo et al. (2007). This method is preferred to measures that make use of the variance of a model as were proposed, for example, by Sobol (1993) or by Saltelli (2002). The variance-based measure requires a high amount of model executions in the sense of Monte Carlo samples that is not affordable in a reasonable time for the model under discussion. On the contrary, the Morris method only requires a small number of executions, which are related to the number of tested parameters and the introduced step size for the parameter modification.
The fundamental idea is the definition of an elementary effect describing the change of the model output y by modifying one of the k parameters p i with a specific step size . Then, the two model outputs y 1 and y 2 are compared to each other. The specific step size = 1∕(l − 1) depends on the number of parameter changes l, called levels. Thus, the parameters p i are usually not in a normalised interval between the normalised lower boundary nlb = 0 and the normalised upper boundary nub = 1 . As a result, has to be transferred to the parameter space defined between the lower and upper boundaries [lb i , ub i ].
The overall method is based on changing one parameter at a time in such a way, that only a minimum of model runs has to be performed and that the parameter space is sampled efficiently. Therefore, a trajectory r within the parameter space is constructed by randomly selecting initial normalised parameter values p i which range from nlb to nub − . Based on this parameter set, one of the parameters p i is then randomly selected and increased by . The whole newly created parameter set itself is then the basis for an increase in another parameter. Further on, this process is repeated until all parameters have changed, cf. Morris (1991) and King and Perera (2013). This results in k + 1 parameter sets. To end up with a sufficient number of elementary effects for a single parameter, r = 10 trajectories are selected resulting in r(k + 1) model evaluations. This is still feasible with regard to the simplified models that are used to identify the crucial parameters. The final quantity of interest in determining the importance of a parameter is achieved by the absolute mean * i of all r absolute elementary effects of a parameter p i yielding cf. Campolongo et al. (2007) and King and Perera (2013). This mean or sensitivity, respectively, allows for a comparison of p i with the other parameters, where the highest values are considered significant, cf. King and Perera (2013). However, this method only allows for an alignment of the relevant parameters, while the decision on what parameters are relevant for the model under discussion has still to be made by the user of the method.
Furthermore, the parameter optimisation strategy is based on maximum likelihood estimations with the goal to determine the model parameters from experimental (110) data, cf. Myung (2003). This is done by relating simulation results y sim to the measurement data y data on the basis of a logarithmic error log = log(y dat ) − log(y sim ) . This error is assumed to be Gaussian distributed resulting in the likelihood function Therein, i is the variance and L the summed likelihood function including both the experimental and the simulated data at point i. In order to optimise the crucial model parameters, the likelihood function L is minimised within the commercial software Matlab 3 and especially within its optimiser fmincon.
For the parameter identification process, two basic experimental procedures are studied, growth by proliferation of cancer cells and atrophy by apoptosis through medical drugs.

Study of the proliferation process
The following section focuses on one of the main cancerrelated processes, the proliferation of cancer cells under sufficient nutrient supply. The corresponding experiments have been carried out by people of the group around Professor Morrison, while the experimental results and data are part of the dissertation thesis of Daniela Stöhr (Stöhr 2018) and a following article (Stöhr et al. 2020).
Despite the numerical model that is able to describe proliferation and atrophy of cancer cells, the model parameters related to these processes have still to be identified and determined. Therefore, the above sensitivity analysis is applied resulting in a model output y composed of the likelihood L of the error between the experimental data and the simulation results. Therewith, the relevant parameters of the proliferation model can be identified and an optimal parameter set can be determined. As a result, the adaptation to the experimental data is obtained.

Proliferation experiment and data
In this study, the human lung cancer cell line NCI-H460 was used that was obtained from the American Type Culture Collection (LGC Standards GmbH, #HTB−177 ). The cells were maintained in RPMI 1640 medium (Life Technologies, Gibco, Karlsruhe, Germany) supplemented with 5% foetal calf serum at 37 • C in a humidified incubator with 5% CO 2 . The experimental procedure was as follows. Firstly, one hundred cancer cells were seeded onto Terasaki multiwell plates for three days and were then transferred onto agarose-coated well plates. Thereon, the cancer cells stuck together and formed a cancer cell spheroid, cf. Figs. 6 and 7. The proliferation of the spheroid was recorded by an inverted digital microscope (EVOS FL Imaging System, Thermo Fisher Scientific, Waltham, MA, USA) and was analysed with the image-editing program Fiji, cf. Schindelin et al. (2012), thus taking the values of the cross-sectional area and the  . 7 Images of the growing tumour cluster at days 3, 7 and 11 after Stöhr (2018) corresponding volume. In a post-processing step, the volume was converted into cell numbers, where the assumption was made that the single cells are spherical occupying a volume of 3.053 × 10 −15 m 3 per cell with no voids in between them. By division of the whole spheroid volume, the amount of cells was obtained. During the experiments, four different sets of twice 11 and once 14, respectively, 15 spheroids have been measured for 13 days. The volume, respectively, the cell numbers, have been summed up and the mean volume m IC vol and the corresponding standard deviation std IC vol have been calculated. N Sp indicates the number of spheroids that have been used for the measurement at the respective day, cf. Table 1. Note that not all spheroids have been measured each day leading to varying spheroid numbers N Sp .

Numerical simulation of the proliferation experiment
Model: In accordance with the experiment, the numerical model aims at determining the proliferation of cancer cells during the time interval of the experiment resulting in a low cancer cell mass as well as volume.
Despite the experimental set-up, where the spheroids are dissolved in a solution, the corresponding numerical simulation treats the spheroid in a simplified 2-d approach. In particular, the spheroid is described here as a certain amount of cancer cells dissolved in an interstitial fluid solution as described in Sect. 2.2.
The basic set-up includes the brain tissue as the solid S together with the interstitial fluid solution I with dissolved nutrients IN and cancer cells IC . Although the blood does not play any role in the proliferation model, the blood is not omitted in this study in order to maintain similarities with the complex overall model and to enable a more reliable transfer between the models. However, B is included in a simplified manner by the use of a constant volume fraction n B = s B n P instead of the angiogenesis formulation based on (67). This leads to a reduced set of governing equations that, in their weak forms, are the overall momentum balance (105), where the gravitational term has been neglected, the volume balance (106) of the interstitial fluid, and the concentration balances (107) of nutrients and cancer cells. This system of coupled equations is solved for the primary variables solid displacement S , blood and liquid pore pressures, p BR and p IR , and nutrient and cancer concentrations, c IN m and c IC m , while the solid volume fraction is obtained from (108).
Finally, the mass interactions (mass production terms) reduce to the proliferation of cancer cells (95) 1 and the nutrient consumption (96).
Parameter identification and optimisation: The goal of this proliferation study is the parameter identification and optimisation of the most sensitive proliferation input parameters. The first step towards this goal is the creation of the model domain. This is handled by the use of the software "Cubit" and results in a 2-d circular finite element (FE) mesh with a diameter of 0.2 m subdivided in 954 Taylor-Hood elements with 2899 nodes, cf. Fig. 8a. The boundary and initial conditions are such that the solid displacements are set to zero along the boundary line, while the external pressures are set to p BR 0 = 1725 Pa and p IR 0 = 258 Pa which, in turn, basically implies a free-flow condition across the boundary for the blood and the interstitial fluid. Furthermore, the nutrient concentration at the boundary is set to c IN m = 1 mol∕m 3 , while the same values of pressures and of nutrient concentration are assumed to be initially distributed in the whole domain. This ensures a sufficient nutrient supply as has been used in the experimental set-up. Furthermore, the cancer cell concentration is restricted to stay within the domain, which is satisfied by the use of a no-flow boundary condition. Initially, the cancer cells occupy the middle of the domain within a radius of 0.0025 m with a concentration of c IC m = 10 −15 mol∕m 3 . This corresponds to a similar amount of cancer cells as in the experiment. Over time, the cancer cells slightly migrate throughout the domain, although the proliferation is mainly concentrated in the middle of the domain, as can be seen from the concentration difference between the beginning and end of the simulation, cf. Fig. 8b, c.
The proliferation model is governed by a lot of input data such as material and model parameters that have a certain sensitivity on the numerical problem under study, cf. Table 2. Therein, the most sensitive parameters are included in the order of decreasing values of the mean * i of their elementary effect. As the values of * i range across 14 orders of magnitude, a horizontal line has been included beneath the most important ones. These are the maximum proliferation rate ⊕ ST, max , the amount N IC of cancer cells per unit mass, the molar cancer cell mass M IC m , the constant K gr indicating half of the maximum of nutrient concentration, and the threshold nutrient concentration c IN m necessary for While all input data of Table 3 are taken from the literature, the most important proliferation data of Table 2 have been optimised by the use of the Morris method, cf. Sect. 5.1. The optimisation was performed for 25 varied starting parameter sets to disturb the maximum-likelihood method in order to find global or at least a variety of different local minima. The best optimised parameter set corresponds to the highest likelihood and is shown in Table 2.
As the molar mass M IC m and the cell amount N IC are reliable quantities, the parameter optimisation basically concerns ⊕ ST, max , K gr and c IN m , cf. Table 4. Note that K SI of Table 2 is the intrinsic permeability of the brain with respect to the interstitial fluid compartment given in m 2 , whereas k I in Table 3 is the corresponding hydraulic conductivity given in m∕s , also compare (49). The same applies to the intrinsic and hydraulic conductivities K SB and k B of the blood.
Comparing the optimisation results with the experimental data as depicted in Fig. 9 reveals, on the one hand, the wide distribution of the single spheroid volume (indicated by purple points) and, on the other hand, the close match of the simulation results (in sky blue) that lay within the standard deviation of the experiment results (purple bars). The results at the individual time steps (in green) illustrate that the time progression follows an exponential function as in the beginning of the cancer cell growth, cf. Laird (1965) and Freyer and Sutherland (1986), as well as in bacterial cell cultures,   cf. Ram et al. (2019), and in yeast, cf. Slavov et al. (2014). A direct comparison of the experimental and the numerical results, cf. Figs. 6 and 9 together with Table 1, exhibits, for example, for day 12, that the optimised curve of Fig. 9 yields a cancer cell volume of 0.19 mm 3 compared to a volume of 0.20 mm 3 given in Table 1.

Study of the apoptosis process
This section focuses on the second main cancer-related process, the cancer cell death triggered by a death-receptor agonist. Thereby, the apoptosis experiments of cancer cell cultures are presented, and the resulting survival-cancer cell ratios are shown. Again, the Morris method is applied in order to find parameter sensitivities and to optimise the most important model parameters. Thus, the significant parameters are adapted to the data and the resulting parameter set as well as a numerical visualisation are given.

Apoptosis experiment and data
The cancer cell survivability in solutions with different concentrations of the therapeutic agent TRAIL, a tumour necrosis factor-related apoptosis-inducing ligand or, more precisely, Db EGFR -scTRAIL, has been experimentally determined using cell staining and flow cytometry techniques, where the material and methods of this experimental set-up are presented in Stöhr (2018) and in Stöhr et al. (2020). In the following, a short summary is given to illustrate the setup. The cancer cells under consideration are again human lung cancer carcinoma cells NCI-H460 that have been grown with 5% foetal calf serum in 2-d cell cultures. The treatment with TRAIL induces apoptosis in the cells by binding to tumour necrosis factor-related apoptosis-inducing ligand receptors (TRAILR1 and TRAILR2) on the cell surface.
In the process of apoptosis, the lipid phosphatidylserine flips from the inside to the outside of the cell membrane, thus allowing staining with Annexin V-EGFP, an enhanced green-fluorescent protein, whose emitted fluorescence can be analysed by the flow cytometry technique. As a result, the percentage of viable cells can be calculated by comparing the number of not stained cells with the overall cell number in a treated cell culture.
During the experiments, the cells have been treated for 24 h with TRAIL resulting in the data shown in Table 5.

Numerical simulation of the apoptosis experiment
Model: With regard to the numerical simulation of the apoptosis process, the overall model is simplified to include only the relevant processes concerning the cell death via TRAIL. The relatively fast cell death during the laboratory experiments allows for major simplifications of the complex overall model. However, despite the experimental set-up of cancer cells bathing within a solution on a plate, the continuum mechanical model is, as for the proliferation experiment before, still related to an initial metastasis ST within brain tissue S . The metastasis is supplied via nutrients IN existing as solutes in the interstitial fluid I . Furthermore, the medication in both the experimental and the numerical setup is applied to the solution that is distributed in the whole interstitial fluid domain. This results in an immediate interaction with the cancer cells or the metastasis, respectively.
To trigger apoptosis, the drug TRAIL is included as a further dissolved component ID of the interstitial fluid. As there is again no angiogenesis as in the proliferation model, the blood is treated likewise. Following this, the simplified set of constituents, solid S with tumour material ST , blood B and interstitial fluid I with dissolved nutrients IN and drugs ID , basically leads to the same reduced set of governing equations as for the proliferation model before. In their weak forms, these are overall momentum balance (105), the volume balance (106) of  the interstitial fluid and the concentration balances (107) of nutrients and drugs, while the solid volume production is computed by the use of (108). This system is again solved for the primary variables, namely the solid displacement S , the blood pressure p BR , the interstitial fluid pressure p IR and the nutrient and drug concentrations, c IN m and c ID m . Parameter study and optimisation: The goal of this apoptosis study is the parameter identification and optimisation of the most sensitive apoptosis input parameters. Here, the same FE mesh is used as in the proliferation study, cf. Fig. 8a or Fig. 10a, respectively. Again, the boundary conditions prescribe no solid displacements along the boundary line together with blood and initial fluid pressures of p BR 0 = 1725 Pa and p IR 0 = 258 Pa and a nutrient concentration of c IN m = 1 mol∕m 3 . The same values are assumed initially distributed in the whole domain, while the drug concentration is similarly applied. Five different scenarios of drug concentrations, namely c ID m ∈ [ 10 −10 , 10 −9 , 10 −8 , 10 −7 , 10 −6 ] given in mol/m 3 , have been investigated at the same nutrient concentration, cf. Fig. 11. When the simulations are initialised, the metastasis is assumed in the centre of the experiment with a diameter of 0.04 m and an initial volume fraction of n ST 0 = 0.1 , cf. the upper part of Fig. 10b. Over a period of 24 h, the drug triggers apoptosis within the cancer cells, thus leading to a cancer cell reduction, cf. Fig. 10c. During the same period, the drug concentrates in the cancerous domain, while the whole domain is shrinking, cf. the lower part of Fig. 10b. As the proliferation model, the apoptosis model is governed by the same default and model parameters, cf. Table 3, with certain sensitivities, cf. Table 6. While all default parameters of Table 3 are taken from the literature, the apoptosis parameters of Table 7 have been optimised as the proliferation parameters before by the use of the elementary effect together with a minimisation of the likelihood function.   [ -]

Fig. 11
The mean viable cancer cell data from the experiment are depicted in purple, while the single viable cell ratios are depicted in dark red. The pink crosses indicate the experiment with the highest values. The simulation results from the optimisation are given in blue including the best result in light green. Additional simulation results corresponding to the best parameter set are given in dark green

Numerical studies
Proceeding from the material and model parameters given in Tables 3 and 8, the following chapter will give insight into the results of numerical experiments on 2-d and 3-d initial-boundary-value problems.
As not all parameters included in Table 8 could be found from the literature or on the basis of experimental data, some of them remained undetermined and had to be marked as chosen. However, those parameters that could be determined by optimisation generally fall into the category of parameters with high values of the absolute mean * i of their elementary effect, cf. Table 6. At the same time, these parameters have a stronger impact on the numerical results compared to those with lower values of * i . For the latter group of parameters, their values had to be chosen, unless they could not be found elsewhere. Although these parameters only have a low impact on the numerical results, this does not mean that they have been chosen arbitrarily. Instead, they have been chosen with respect to sound numerical results. For example, the parameters B , angio , A , B and c IV m together with B, angio after Eq (70) govern the relation of the sigmoid function (67). These values are chosen for a reasonable steepness of the functions s B (p dif ) , p dif (s B ) and BR B (s B ).
In addition to the above and in contrast to the values given in Table 3, the Lamé constants have been changed towards the values of Table 9. While the parameters of Table 3 represent white matter parameters ) that correspond to the experimental set-up, the values of Table 9 are based on an optimisation of the whole brain, white and gray matter, cf. Soza et al. (2004), and additionally emphasise values that have been found by experiments on pig brains.
The following studies proceed again from the full set of governing equations (105)-(108) that are further solved by the use of the FE solver PANDAS. This procedure guarantees a fully coupled numerical analysis of all features including growth and atrophy of cancer-cell clusters as well as solid deformation and stress, flow of blood and the interstitial-fluid mixture together with diffusion of fluid species.

Two-dimensional study of the overall model
Based on the above statements, a 2-d computation of growth and atrophy of a single lung cancer metastasis within the brain tissue is simulated. The goal of the following study is to present the capability of the introduced model including cancer cell proliferation, angiogenesis, necrosis and apoptosis. While the numerical examples of the preceding chapter, that have been needed to identify material parameters, exhibit boundary-value problems with rotational symmetry, the actual problem is fully two-dimensional.

Set-up
The model includes all introduced primary variables, namely the solid displacement S , the blood and interstitial fluid pressures, p BR and p IR , as well as the nutrient, cancer cell, VEGF and drug concentrations, c IN m , c IC m , c IV m and c ID m . As before, the study is carried out on a circular domain with a radius of 0.2 m representing a geometrically simplified model of a brain bathing in a rigid skull with an initial cancer cell cluster in the middle, cf. Fig. 12. Concerning the initial and boundary conditions, the solid skeleton is fixed along the boundary with ̄ S ≡ . Constant fluid pressures p BR 0 = 1725 Pa and p IR 0 = 258 Pa are prescribed at the boundary and initially in the entire domain together with n B 0 = 0.05 and n I 0 = 0.20 . Moreover, all fluid fluxes across the boundary are set to zero with the exception of the interstitial fluidsolvent flux ̄ IL ⋅ . In addition, the injected interstitial fluid flux ̄ I ⋅ including the drug-species at the yellow point in the middle of the left side, cf. Fig. 12, can vary over time, such that a therapeutic infusion can be applied.
The infusion is assigned to last for 8 h including a 2-h period of linear infusion rate increase and another 2-h period of linear infusion rate decrease. Altogether, two infusion scenarios are investigated. In the first case, the maximum medication rate reaches 1.5 × 10 −10 mol∕(m 2 s) at the infusion site resulting from a therapeutic agent concentration of 1.5 × 10 −3 mol∕m 3 at an interstitial fluid velocity of 10 −7 m∕s . In a second case, the maximum medication rate is reduced to 1.5 × 10 −11 mol∕(m 2 s) resulting from 1.5 × 10 −4 mol∕m 3 at the same interstitial fluid velocity of 10 −7 m∕s . As the medical drug is injected with the interstitial fluid, the dissemination of the drug is governed not only by diffusion but also by the fluid flow, thus resulting in a convection-enhanced delivery (CED) of the medical drugs. Furthermore, note in passing that in contrast to previous studies, such as Linninger et al. (2008) or Ehlers and Wagner (2015), the infusion site is in this example not located at the end of an infusion needle intruded into the brain but at the outer boundary. As a result, the motion of the interstitial fluid with included medical drug is governed by the fluid's influx velocity of 10 −7 m∕s . In case of an approximately vanishing solid velocity, the influx velocity equals the seepage velocity, such that the corresponding filter velocity results in 2 × 10 −8 m∕s at an interstitial fluid volume fraction of n I = 0.2.

Results
For the initial-boundary-value problem described above, Fig. 13 exhibits the process of proliferation and related nutrient consumption as well as the response of a growing tumour to an early treatment. The latter, of course, is only possible if the tumour has been discovered at this early state before angiogenesis has started. The results presented in Figs  Proliferation: Once a metastasis-based tumour cell cluster has settled in the brain, a volume increase of the tumour takes place under fully nutrient-supplied conditions. This can be seen from Fig. 13, where the tumour volume fraction n ST , the tumour mass m ST , the nutrient concentration c IN m and the drug concentration c ID m are plotted versus time. It is furthermore seen from Fig. 13a-c that the tumour increases exponentially beyond the micrometastatic switch at day 58, while the nutrient concentration strongly decreases as a result of nutrient consumption of the growing tumour. This behaviour stops at day 77, when the tumour has been detected and a medical treatment has been initiated, cf. Fig. 13d.
Apoptosis: Once the medication has started, the dark and light blue lines follow different paths describing the situation at Point A of the tumour site. While the dark blue line corresponds to a maximum drug infusion rate of 1.5 × 10 −10 mol∕(m 2 s) , the light blue line proceeds from a maximum value of only 1.5 × 10 −11 mol∕(m 2 s) . Following this, it is seen from Fig. 13a that, depending on the medication rate, the dark blue line exhibits a much steeper descent of the tumour volume fraction than the light blue one. Furthermore, when the medication stops at approximately day 80, the tumour recreates much faster when the lower medication has been applied instead of the higher one. These results also correspond to the nutrient concentration at the tumour site that decreases when the tumour grows and recreates when the tumour growth ends, however, following different paths depending on the chosen drug infusion rate.
Corresponding to the course of the tumour volume fraction n ST , one observes the growth of the aggregated tumour mass m ST , cf. Fig. 13c, that reaches 40 mg when the medication starts. It is furthermore seen from Fig. 13d, how different the two medications represented by the dark and light blue lines are. While the dark blue line describes a drug concentration of the infusion of 1.5 × 10 −3 mol∕m 3 , the light blue line only results from a concentration of only 1.5 × 10 −4 mol∕m 3 . As a result of the CED process, the drug concentration at Point A at the tumour site, cf. Fig. 13d, reaches approximately 1.2 × 10 −7 mol∕m 3 , in the first case, and approximately 1.2 × 10 −8 mol∕m 3 , in the second case.
Angiogenesis and necrosis: When the medication has stopped at day 80, the tumour starts regrowing, cf. Fig. 13a. Note that this regrowth cannot be observed in Fig. 14c as a result of the scale difference between these two figures. However, it is seen from Fig. 13a that, depending on the medication rate, the dark blue line exhibits a much steeper descent of the tumour volume fraction than the light blue one. Furthermore, after the medication, the tumour recreates much faster when the lower medication has been applied instead of the higher one. These results also correspond to the nutrient concentration at the tumour site that decreases when the tumour grows and recreates when the tumour growth ends, however, following different paths depending on the chosen drug-infusion rate. On the other hand, the nutrient consumption belonging to the dark blue line is so high that the dark blue line touches the purple line, while the light blue line enters the critical region at the umber line. However, the tumour succeeds in both cases to continue growing.
Stimulated by the heavy nutrient consumption, VEGF is released at day 84 and blood vessels start to grow towards the tumour, which also results in an increase of the nutrient concentration, cf. Fig. 14b, d. After day 120, the nutrient concentration decreases again and the tumour would further on undergo necrosis unless VEGF is released again. Figure 15 exhibits the volumetric solid deformations expressed by the Jacobian determinants J S of the total deformation, J Sg of the deformation induced by growth (including apoptosis), and J Sm of the purely mechanical deformation, cf. (73) and (74). Solid lines correspond to the higher medication of 1.5 × 10 −3 mol∕m 3 and dashed lines to the lower medication of 1.5 × 10 −4 mol∕m 3 with J S in blue, J Sg in green and J Sm in purple. As the Jacobians stay at 1.00 in the undeformed state up to approximately day 100, it is concluded that measurable deformations only start after the angiogenesis. Figure 15a displays Jacobians with values either greater or less than one resulting in volumetric extensions or compressions. As the stresses are obtained from J Sm = J S ∕J Sg in purple, the growing tumour leads to compressive stresses after angiogenesis, cf. Fig. 15a. Figure 15b that, as well as Fig. 15d represents its values at day 125 along a horizontal line through Point A, moreover exhibits that the highest values of the Jacobian determinants are not at the centre of the domain at 0.1 m but shifted to the left and right. On the other hand, the centre of the metastasis between the two peaks exhibits necrosis leading to a reduction of the volume fraction, whereas the boundary could be supplied with nutrients and proliferates as a result of the angiogenesis process. Figure 15c shows the von-Mises stresses with higher values for the lower and lower value for the higher medication. Furthermore, Fig. 15d presents the same shifting behaviour for the von Mises stresses as for the Jacobians with the highest value of approximately 277 Pa. Please note that small numerical instabilities occurred that can be seen in these figures due to the sharp switch based on the Heaviside function included in the apoptosis model (94) and (99) and the von-Mises stress, cf. Fig. 15, it is seen from Fig. 16 that the effective volumetric stress T S ef f , vol as well as the partial pore pressures s B p BR and s I p IR only reacts to the tumour growth after angiogenesis initialised around day 83. With increasing growth, the brain tissue is pressed aside and gets under pressure. Once this has happened, it is seen from Fig. 16a that the reaction to tumour growth at Point A strongly depends on the medication with the stronger reaction depending on the lower drug medication (dashed lines). Figure 16b displays the volumetric stress and the pressures along a horizontal line through Point A. Here, it is seen as a result of the injection site at the left of the domain that at day 125 the tumour centre has been shifted to the right with an amount depending on the medication. Note that necrosis occurs in both cases, but is not as dominant as apoptosis due to the medication. However, as the tumour regrows after the medication has stopped, the effect on stress and pressures is larger in case of the lower medication. In addition to the individual graphs included in Fig. 16a, b, it is worth to compare T S ef f , vol to the effective pore pressure p FR , thus comparing the two parts of T vol to each other. This can be done by taking the ratio of T S ef f , vol over the negative effective pore pressure yielding As p FR = s B p BR + s I p IR is nearly constant, p FR in (114) has been substituted by its initial value p FR 0 with the result that Figure 16c exhibits the course of r vol at day 125 along a horizontal line through Point A with peak values of 10.80% (high medication) and 32.08% (low medication). Finally, it is again seen that the tumour centre shifts to the right depending on drug concentration of the infusion. To close up, Fig. 17a shows the metastasis volume fraction before the medical treatment has been started with the tumour site in the centre of the domain, while Fig. 17c exhibits the same tumour, however, once the treatment is over and the drug concentration c ID m has decreased below the threshold level of c ID m = 10 −11 mol∕m 3 . As before and as result of the infusion from the left side, it is seen that the centre of the reduced metastasis is now more to the right than before. Completing the 2-d example, Fig. 17b displays a state in between the situations at (a) and (c), where the therapeutic agent concentration is given together with the arrows of the seepage velocity IL of the interstitial fluid solvent with colour-coded normalised values. Note again that, as a result of the little masses of the fluid solutes compared to the solvent mass, IL ≈ I .  Soza et al. (2004) 6.2 Three-dimensional simulation of a tumour treatment

Set-up
The present example, cf. Fig. 18, is based on a 3-d domain with a diameter of 0.2 m and a red mark in the centre. In a realistic scenario, a brain tumour metastasis is usually detected by sufficiently resolved medical imaging techniques, once it has reached a certain size. In the present example, the tumour is assumed to have been found as the irregular white region around the centre of the domain, such that a surgical intervention can be planned. The initial location and shape of the cancer cells or of the metastatic volume fraction, respectively, are based on a MRT scan of a meningioma which is used here in the sense of a proxy with data originating from the 3DSlicer Registration Case Library found on the MIDAS webpage (https:// www. insig ht-journ al. org/ midas) in the section National Alliance for Medical Image Computing (NAMIC). Apart of the spherical domain consisting of 5104 spatial Taylor-Hood elements, the problem is completely non-symmetric, meaning that a fully 3-d problem had to be investigated. The initial and boundary conditions are chosen as the same that have been applied in the 2-d problem before. The planned surgery has been defined as a CED-based medical drug infusion through a little whole in the skull, such that an infusion needle with inner radius of 3.6 mm could be inserted at a length of 50 mm. Different infusion scenarios, cf. Fig. 20, have been studied in order to find the most appropriate surgery result.

Results
During neurosurgical drug injection of the metastasis, cf. Fig. 19, four infusions (Phases I − IV ) with a drug concentration of c ID m = 2.5 × 10 −4 mol∕m 3 are given at an interval of five days. Numerically, Phase I displays the treatment 25 s after the computation has been initiated with the tumour region around Point A in black. Obviously, the highest drug concentration is here exactly at the outlet of the infusion needle (Area I) surrounded by red and purple areas with threshold concentrations of 5 × 10 −8 mol∕m 3 and 10 −11 mol∕m 3 . Once the therapeutic drug reaches the tumour site, the degradation of tumour material starts displayed by a colour change of the tumour volume fraction from black to white. After the second treatment five days after the first infusion (Phase II), the whole domain is at least in purple with a drug concentration of more Spherical domain and mesh of the 3-d initial-boundary-value problem (diameter 0.2 m) with boundary conditions as in the 2-d problem before, tumour site seen as black area with a red mark (Point A) in the centre, and medical treatment through the infusion needle (inner radius 3.6 mm, insertion length 50 mm) close to the top than 10 −11 mol∕ m 3 , while the red area with c ID m ≥ 5 × 10 −8 mol∕m 3 has reached the tumour. As a result, the degraded part of the tumour is rapidly growing and continues growing after the third treatment 10 days after the initiation of the treatment (Phase III). This behaviour continues after the forth treatment 15 days after the first treatment at Phase IV 1 , such that the tumour has nearly vanished half a day after the forth infusion (Phase IV 2 ). Only a small region at the bottom of the tumour remains. However, exactly this is the problem. Without further medication, if medically meaningful, the tumour regrows and reaches at the end of the computation at t end ≈ 30.7 days later the displayed condition. Figure 20 exhibits numerical results obtained for different infusion scenarios. In particular, Fig. 20a, b exhibits the degradation and regrow of a brain metastasis under different TRAIL applications. Based on a single or a twice applied drug infusion, it is obvious that a higher drug concentration results in a stronger degradation of the tumour, although a regrow could not be avoided. In addition to this, Fig. 20c, d, presents multiply applied drug infusions with different TRAIL concentrations. The numerical studies presented in this chapter provide an option to plan a surgical treatment in advance, such that an optimal treatment with the most suited infusion protocol can be offered to the patient.

Concluding remarks
In the present article, the computation of metastatic lung cancer cell proliferation and atrophy in brain tissue has been investigated based on experimental data that have been obtained by people around Professor Morrison at the University of Stuttgart. Based on these data, it was possible to calibrate a continuum mechanical model that has been elaborated on the basis of the Theory of Porous Media (TPM). Only by the use of the TPM, it was feasible to include a deforming solid skeleton, the porous solid brain structure, an interstitial fluid and blood as mutually acting constituents of the overall model. As living tissues need nutrients and further ingredients for their basal nutrition and growth, the interstitial fluid has been treated as a real mixture of a fluid solvent and various solutes or species, respectively. When cancer comes into play, a clustered cell structure nests in the skeleton tissue, here the brain, and starts to grow without accepting apoptosis signals. The fact of growth, in particular proliferation and atrophy, has been included in the description by the introduction of a multiplicative decomposition of the solid deformation gradient in a mechanical and a growth-dependent part, on the one hand, and by the addition of a mass production term to the mass balance of the solid skeleton that interacts with the interstitial fluid, on the other hand. Based on the integration of production terms in the mass balance equations of solid and interstitial fluid, it could be shown that this procedure leads to the definition of a growth-dependent intermediate configuration that is different from the usual definitions in elasto-plasticity or in solid-based thermomechanics.
The constitutive equations of the overall model that had to be formulated for the complex description of solid deformation mutually coupled with flow processes of the interstitial fluid and the blood have been proven to fulfil the thermodynamical restrictions resulting from the entropy inequality of the multiphasic aggregate. In conclusion, the full set of constitutive equations including the fluid species that trigger cancer proliferation or apoptosis, the latter in the course of a medical treatment, is convenient for the computational description of the full and complex model under consideration. Therewith, a model has been designed that can be used for pre-clinical studies to further validate and optimise model performance and to ultimately improve treatments and treatment outcomes to optimise the application of medical drugs, for example, in the skull. (120) div � = (det ) � det .
When isotropic growth comes into play, as it is the case in the present article as a result of tumour growth, the above procedure can also be applied to growing material yielding As growing material is governed by the general mass balance, cf. for example (129), this gives rise to the a priori constitutive equation such that where det S = det Sm det Sg has been used, thus furthermore substituting (128) in case of materially incompressible S . From (135), it is seen that (det Sm ) −1 plays the role of a push-forward transformation of n S 0 belonging to the reference configuration S 0 towards n S belonging to the current configuration S . From (135), we additionally have where both, det S and det Sg , are pull-back transformations of volume fractions pulling n S from S and n S 0 from S 0 towards n S g of the same configuration S g . As a result, Therein, det S , det Sm and det Sg do not only define pullback transformations of volume fractions but also push-forward transformations of volume elements yielding where the corresponding deformation gradients S , Sg and Sm transform line elements via By these means, one obtains the order of configurations as in Fig. 21 differing from the standard literature, cf. Sect. 3.3.
n S g = n S det S = n S 0 det Sg , while n S 0 = n S det Sm .
(138) dv = det S dv g , dv 0 = det Sg dv g , dv = det Sm dv 0 , (c) Nonlinear elasticity of growing brain tissue As in thermoelasticity, the purely mechanical deformationversus-stress behaviour is governed by Sm also defining the right Cauchy-Green tensor Sm , cf. (77) 1 , while the growth behaviour depends on (133). Based on this equation, the growth-depending Cauchy-Green strain yields As S , Sm and Sg are coupled to each other through the multiplicative decomposition (133), only two of them are independent, while the third is a function of these two. In terms of the exploitation of the entropy inequality of the whole system, cf. (33), it is convenient to proceed with S = T S S and Sg , where S is directly coupled to the solid motion trough the solid displacement S and the deformation gradient S , while Sg is a function of growth through n S or ̂S , respectively, cf. (134).
From the entropy inequality, one has (41) reading where, by means of the principle of determinism, S was assumed as a function of S . Based on the principle of frame indifference, the dependence of S on S must be found in such a way that S depends on S through S . This results in (78). To obtain (78)