Sparse optimal control of a phase field tumour model with mechanical effects

In this paper, we study an optimal control problem for a macroscopic mechanical tumour model based on the phase field approach. The model couples a Cahn--Hilliard type equation to a system of linear elasticity and a reaction-diffusion equation for a nutrient concentration. By taking advantage of previous analytical well-posedness results established by the authors, we seek optimal controls in the form of a boundary nutrient supply, as well as concentrations of cytotoxic and antiangiogenic drugs that minimise a cost functional involving mechanical stresses. Special attention is given to sparsity effects, where with the inclusion of convex non-differentiable regularisation terms to the cost functional, we can infer from the first-order optimality conditions that the optimal drug concentrations can vanish on certain time intervals.


Introduction
Mechanical stresses play a significant role in both enhancing and inhibiting the growth of tumours. The unregulated proliferation of tumour cells displaces nearby normal tissues and in turn these tissues exert externally applied stress to resist tumour expansion. In various experimental studies (see [9,25,26,44] and the references cited therein) high compressive stress has the effect of suppressing proliferation and can induce apoptosis (natural cell death) in tumour cells. However, in the case where the mechanical loads are not uniform, tumours can adapt by growing in directions of least stress. Moreover, deformations of the microenvironment brought about by these mechanical loads can alter the structure of nearby blood and lymphatic vessels, which are responsible for supplying the region with crucial nutrients, oxygen, therapeutic drugs, as well as drainage of excessive interstitial fluids containing waste products. The gradual reduction in blood flow turns the stressed region more hypoxic and more acidic, compounding with the reduction in nutrients levels further accelerates the invasive and metastatic potentials of the tumours cells. On the other hand this also impairs the effectiveness of immune cells or therapeutic agents as they are not able to reach certain tumour regions in sufficient quantities. With the use of mathematical modelling [29,43], treatments aimed at alleviating stress seem to be a promising avenue that warrant further investigations and could be used in coordination with other anti-cancer therapies.
Recent progress in mathematical oncology have shown promising results in forecasting tumour growth and predictive simulations of treatments [1,2,32,33,34,35]. Most models employ a continuum description involving partial differential equations to capture a multitude of biological and chemical mechanisms. Among those we focus on the subclass of phase field tumour models [13,23,38,45], where the corresponding numerical simulations (see, e.g., [14,15,17,45,46]) are able to replicate commonly observed morphologies exhibited by tumours and their vasculatures.
While there has been a surge of activity in the subsequent mathematical modelling and analysis of phase field tumour models, see [1,11,14,15,23,38,46] and the references cited therein, there seems to comparatively fewer focus on mechanical interactions in tumour growth within the subclass of phase field models, aside from recent contributions [16,22,32,33], see also [3,5,19,20] for results concerning the related Cahn-Larché system. In light of the significance of mechanical stress, for our study, we consider a simplification of the phase field model that was proposed and studied in the authors' previous work [22]. Consider a bounded domain Ω ⊂ R d , d ∈ {2, 3}, with boundary Γ := ∂Ω that is either of C 1,1 -regularity or is convex and is partitioned into two subregions Γ D and Γ N . For an arbitrary T > 0 (which can be interpreted as the length of the medical treatment), the following model posed in the space-time cylinder Ω × (0, T ) describes the evolution of a cellular mixture containing tumour and non-tumour cells subject to various mechanisms involving a chemical species acting as nutrient and mechanical stresses: We refer the reader to [22,32,33] for more background on the model and related topics, while briefly describe the main components. In the above, the variable ϕ denotes a phase field parameter that serves to distinguish between the two different types of cellular material in the mixture, with tumour cells occupying the region {ϕ = 1} and non-tumour cells in the region {ϕ = −1}. The subsystem (1.1a)-(1.1b) constitutes a Cahn-Hilliard type equation, where µ is the associated chemical potential. Coupled to this is a reactiondiffusion equation (1.1d) for a nutrient σ, as well as a quasistatic linear elasticity system (1.1e) with displacement u and symmetric strain tensor E(u) := 1 2 (∇u + (∇u) ⊤ ). We mention that there are certain cases where the nutrient evolves quasistatically, which is covered by the case β = 0. The terms W ,ϕ (ϕ, E(u)) and W ,E (ϕ, E(u)) are partial derivatives of the elastic energy W(ϕ, E(u)) with respect to its first and second arguments, respectively, and for this work we consider the choice where C is a constant, symmetric and positive definite elasticity tensor satisfying the usual symmetry conditions, and the phase-dependent stress-free strainĒ(ϕ) under Vegard's law is given by the linear ansatzĒ(ϕ) =Ē + ϕE * with constant symmetric second order tensors E and E * . Furthermore, in (1.1b) the directed movement of cells by chemotaxis is captured by the term −χσ, with χ ≥ 0 playing the role of chemotactic sensitivity [23], while the term Ψ ′ (ϕ) is the derivative of a double-well potential Ψ(ϕ) with equal minima at ϕ = ±1. In our setting this term plays the role of cellular adhesion that leads to the development of regions of tumour and non-tumour cells well-separated by interfacial layers described by the set {−1 < ϕ < 1}. For boundary conditions we subdivide the boundary Γ into the partition Both portions are assumed to be relatively open and to have positive Hausdorff measures, and on the portion Γ D , representing a rigid structure of the tumour environment such as bone, the displacement u is set to be zero, and on the complement portion Γ N , the normal component of the stress tensor W ,E is equal to some given load g provided by a fixed source. Meanwhile, (1.1h) highlights that the cellular diffusive flux ∂ n µ is zero across the boundary, and for κ > 0 the nutrient flux ∂ n σ is proportional to the difference between a nutrient source σ B from nearby capillaries and the nutrient level at the boundary. The case of a zero nutrient diffusive flux is covered by the choice κ = 0. Lastly, the source term U (ϕ, σ, E(u)) in (1.1a) captures cellular growth that can be influenced by nutrient concentration and mechanical stress. The example we will use is where λ p ≥ 0, λ a ≥ 0 are constant proliferation and apoptosis (cell death) rates, and f , g and k are Lipschitz and bounded functions. For instance, we can model the proliferation and apoptosis of tumour cells only by prescribing the conditions f (1) = k(1) = 1, f (−1) = k(−1) = 0, see, e.g., [23]. Meanwhile, to account for the effect of reduced proliferation due to increase in mechanical stress [4,9,25,44], we may consider as a motivating example the function g : R d×d → R defined as where |A| is the Frobenius norm of the matrix A, so that as the magnitude of the stress W ,E (ϕ, E(u)) increases, the effects of the proliferation term λ p σf (ϕ)g(W ,E (ϕ, E(u))) become less significant. This is different to the choice considered in [22] as the derivation of optimal conditions in our present contribution requires a differentiable g. What is not present in the previous work [22] is the coefficient m(t), and when paired with k(ϕ) we use the product m(t)k(ϕ) to model a cytotoxic drug-induced decrease in tumour proliferation. A motivating example for m(t) from [11] is with drug dosage d c , drug delivery times T i for i = 1, . . . , n, where n is the number of chemotherapy cycles, τ denoting the mean lifetime of the drug and H is the Heaviside function. After the i'th infusion, the effect of the drug decreases exponentially until the next infusion at time T i+1 . For drugs with sufficiently short mean lifetime τ , or with large enough infusion gap T i − T i−1 , there are certain time intervals where the coefficient m is close to zero. Similarly, the source term S(ϕ, σ) in (1.1d) accounts for nutrient consumption and transport to and from external capillaries. The example we will use is of the form with constant consumption rate λ c ≥ 0, capillary supply rate B ≥ 0, capillary nutrient concentration σ c and a Lipschitz, bounded function h. For instance, we can model nutrient consumption only by the tumour cells by prescribing the conditions h(1) = 1 and h(−1) = 0. A new element absent from [22] is the coefficient s(t), which models the reduction in nutrient supply caused by antiangiogenic therapy, and in [11] a similar form to (1.3) is proposed for s(t), meaning that under suitable conditions, the coefficient s(t) take values close to zero for certain time intervals.
It is common to prescribe cytotoxic drugs in chemotherapy that serve to disrupt the cellular division process and promote apoptosis, but tumours can overcome these effects by developing drug resistance or by generating new vasculatures through angiogenesis to obtain nutrients that compensate any loss of mass from chemotherapy. Therefore, in certain situations, it is of interest to combine two or more different therapies so that their joint effect can account for more mechanisms that allows tumours to avoid complete elimination, and have an overall larger positive impact on the treatment than the individual monotherapies. Unfortunately, the results of various experimental and clinical studies (see [36] and the references cited therein) have not produced clear guidelines on how to proceed with combined therapies, in part due to the multitude of drugs presently available and patient-specific interactions of multiple drugs. Hence, mathematicians and physicians have turn towards the framework of optimal control to infer protocols, dosages and timings that maximise tumour reduction and minimise harmful side-effects [27,30,31,37,39]. To contribute to this effort, we study an optimal control problem with the model (1.1) as the state system, and as controls we work with the boundary nutrient supply w 1 = σ B , the cytotoxic coefficient w 2 = m(t) and the antiangiogentic coefficient w 3 = s(t). The cost functional we consider is It is composed of the standard tracking-type with weights α Q , α Ω ≥ 0 and target functions ϕ Q : Q → R and ϕ Ω : Ω → R, and L 2 -regularisations for the optimal controls w 1 = σ B , w 2 = m(t) and w 3 = s(t) with corresponding weights γ 1 , γ 2 , γ 3 ≥ 0. Let us stress that the controls w 2 and w 3 are solely functions of time and are spatially constant. Compared to previous works on the optimal control with phase field tumour models, we have the presence of a term involving the the square of the stress W ,E (ϕ, E(u)) weighted by a nonnegative coefficient n(x, ϕ) and constant α E ≥ 0. Due to the role of mechanical stresses on enhancing tumour growth, we are motivated to minimise stress accumulating in a certain region of the domain, such as important organs (by taking n(x, ϕ) = χ D (x) for a subregion D ⊂ Ω where χ D is the characteristic function of the set D), or in certain parts of the tumour microenvironment whose location can be encoded with the help of the phase field variable ϕ. One example is a function n(x, ϕ) = max(0, min(1, 1 2 (1 − ϕ))), so that n is non-zero in the non-tumour region {ϕ = −1} and is zero in the tumour region {ϕ = 1}.
Moreover, we prescribe L 1 -regularisations of the drug concentrations w 2 and w 3 with weights γ 4 , γ 5 ≥ 0 to the cost functional (1.4), with the aim of using the combination of both L 2 and L 1 -regularisations to show sparsity, see Theorem 7 below for the precise formulation. A first work on sparse controls with phase field tumour models is [41], where directional sparsity [24] of the controls, i.e., sparsity w.r.t. space or w.r.t. time, is shown. Our reasoning for such considerations is in part motivated by the common practice that chemotherapies should be administrated to the patient only in very short periods of time to avoid adverse side-effects. In the simulations performed in [12], where an optimal control problem of a similar nature is studied with only L 2 -regularisation terms in the cost functional, the optimal cytotoxic drug concentration is positive over the treatment period. In practical applications this translates to prolonged exposure and subsequent accumulation of the drugs in the body, potentially invoking damaging side-effects and may even entail a premature abortion of the medical treatment.
The goal of this paper is to study the optimal control problem (1.4) subjected to the state system (1.1). Building on the well-posedness results established in [22], we prove the existence of a minimiser and derive first-order optimality conditions. Our main result is sparsity of the optimal drug concentrations, brought about by the convex non-differentiable L 1 -terms in (1.4). Compared to [41], our analysis includes the elasticity interactions in (1.1b) and in (1.4), covering both cases of β > 0 and β = 0 in (1.1d) in a uniform manner, as well as different sparsity conditions for non-negative drug concentrations m(t) and s(t).
We comment that tracking terms involving the nutrient concentration σ, such as can also be inserted into the cost functional. Other terms of interest include the total tumour volume at time T given by the spatial integral of 1 2 (1 + ϕ(T )), and thanks to the well-posedness result for (1.1) (see Theorem 1 below) we can consider other parameters as control variables, for instance the capillary nutrient concentration σ c , the boundary load g, the initial data ϕ 0 , σ 0 and the coefficients χ, λ p , λ a , λ c in (1.1) in the context of parameter estimation [18,28], and even the magnitude of the treatment time T [8,21,40]. One can also consider spatially varying drug concentrations m(t, x) and s(t, x) as in [12,41], and the corresponding analysis to adapt to these elements would only require minor and straightforward modifications.
The rest of the paper is organised as follows: We recall previous results in Section 2, and the existence of a minimiser to the optimal control problem is shown in Section 3. Section 4 is devoted to the derivation of first-order optimality conditions, and in Section 5 we discuss the sparsity of controls.
2 Mathematical setting and previous results

Notation and useful preliminaries
The standard Lebesgue and Sobolev spaces over Ω are denoted by L p := L p (Ω) and W k,p := W k,p (Ω) for any p ∈ [1, ∞] and k > 0, with corresponding norms · L p and · W k,p . In the case p = 2, these become Hilbert spaces and we use the notation H k := H k (Ω) = W k,2 (Ω) and the norm · H k . For any Banach space Z, we denote its dual by Z ′ , and the corresponding duality pairing by ·, · Z . When Z = H 1 (Ω), we use the notation ·, · = ·, · H 1 . The L 2 (Ω)-inner product is denoted by (·, ·), while the L 2 (Γ) and L 2 (Γ N )-inner products are denoted by (·, ·) Γ and (·, ·) Γ N , respectively. We define the Sobolev space H 2 n (Ω) as the set {f ∈ H 2 (Ω) : ∂ n f = 0 on Γ}, and for the displacement u, we introduce the following function space: where by [10, Thm. 6.15-4, pp. 409-410], a Korn-type inequality is valid in X(Ω), i.e., there exists a constant C K > 0 such that

Assumptions and previous results
In this work we make the following assumptions regarding parameters and functions in the model: (A1) Let g ∈ L 2 (Γ N ) d and σ B ∈ L ∞ (Σ) be given, while β, B, κ, χ, λ a , λ p , λ c , σ c are nonnegative constants such that at least one of {B, κ} is non-zero if β = 0. Moreover,Ē and E * are constant symmetric second order tensors while C is a constant symmetric, positive definite fourth order tensor satisfying for all symmetric second order tensors E ∈ R d×d sym with a positive constant c 0 .
(A2) The potential Ψ = Ψ 1 + Ψ 2 is a non-negative function, Ψ i ∈ C 3 (R) for i = 1, 2, with a convex non-negative function Ψ 1 such that for all r, z ∈ R, for some positive constant C.
Lipschitz constants that shall be denoted by a common symbol L > 0. Furthermore, we assume h is non-negative.
It is worth noting that the conditions expressed in (A2) are fulfilled by the classical quartic potential Ψ(r) = 1 4 (r 2 − 1) 2 . For the motivating example (1.2), for any A ∈ R d×d , we use the notation g ′ (A) to denote the tensor derivative of g, i.e., g ′ (A) is a second order tensor with On the other hand, we use the notation g ′′ (A) to denote the Hessian of g, which is a fourth order tensor defined as Hence, it is easy to see that for any In particular, we can infer that |g ′ (W ,E (ϕ, E(u)))| and |g ′′ (W ,E (ϕ, E(u))| are bounded a.e. in Q thanks to (A3). For the rest of the paper, the parameters β, g, B, κ, χ, λ a , λ p , λ c , σ c , C,Ē, E * , as well as initial data ϕ 0 and σ 0 are kept fixed. We then introduce the notation and the set of admissible controls U ad = U The admissible set of controls U ad is a non-empty, closed and convex subset of U := L 2 (Σ) × L 2 (0, T ) × L 2 (0, T ), and we can find a positive constant R such that The following result concerns the well-posedness of the model (1.1).

The optimal control problem
In this section we show that there exists at least one minimiser to the optimal control problem minimising the cost functional (1.4) with state system given by (1.1). By Theorem 1, we can define the control-to-state operator S which assigns every admissible control w = (w 1 , w 2 , w 3 ) = (σ B , m, s) the corresponding unique solution (ϕ, µ, σ, u) to (1.1). Namely, we have where the solution space Y β is defined, according to Theorem 1, as Denoting by S 1 (w) = ϕ the first component and by S 4 (w) = u the fourth component, we can define the reduced cost functional as Since the proof is somewhat standard we omit the details and sketch the main points. The non-negativity of J implies the infimum inf U ad J exists and allows us to find a minimising sequence {w n = (w 1,n , w 2,n , w 3,n )} n∈N ⊂ U ad such that J (w n ) → inf U ad J as n → ∞. Denoting the corresponding solution as (ϕ n , µ n , σ n , u n ) = S(w n ) ∈ Y β , we infer by the bound (2.4) that {(ϕ n , µ n , σ n , u n )} n∈N is uniformly bounded in Y β . Hence, along a non-relabelled subsequence there exists a limit triplet w = (w 1 , w 2 , w 3 ) ∈ U ad such that, as n → ∞, The Aubin-Lions compactness theorem then yields the strong convergence of ϕ n to ϕ in C 0 ([0, T ]; L 2 (Ω)), allowing us to pass to the limit in the tracking terms of J and provides strong convergence n(ϕ n )η → n(ϕ)η for all η ∈ L 2 (Q). Together with the weak convergence of E(u n ) to E(u) in L 2 (Q), we arrive at the weak convergence n(ϕ n )W ,E (ϕ n , E(u n ))→ n(ϕ)W ,E (ϕ, E(u)) weakly in L 2 (Q).
The proof of Theorem 3 proceeds in four steps, which is covered in the following four subsections.
Proof. We proceed with formal estimates that can be justified rigorously with a Galerkin approximation. In the following the positive constants denoted by the symbol C will be independent of the Galerkin parameter, as well as h 1 , h 2 and h 3 and might change from line to line. Besides, let us remark that since w * is fixed, the corresponding state (ϕ, µ, σ, u)= S(w * ) enjoys the bound (2.4). Let us mention that uniqueness follows from existence thanks to the linearity of the system (4.3). We test (4.4a) with η and Kξ, (4.4b) with −ξ t and η, (4.4c) with Rψ and (4.4d) with v t for positive constants K, R to be determined later. After summing and rearranging we get 1 2 and we have used the identity derived with the help of the symmetry of C and E * . Furthermore, by the positive definiteness of C and Young's inequality, Next, recalling that w * 2 and w * 3 are constant in space, from the definition of U lin and S lin , using the boundedness of w * 2 , w * 3 , f , g, h, k and their derivatives, as well as the boundedness of σ, and of |g ′ (W ,E (ϕ, E(u)))|, it is easy to see that while where we also use that the first two terms on the right-hand side are non-positive. Also, using (A2), the continuous inclusion V ⊂ L 6 (Ω), and the fact that ϕ ∈ L ∞ (0, T ; H 1 (Ω)), Hence, the right-hand side of (4.5) can be estimated as follows: with a positive constant c that is independent of R. To handle the term involving ξ t , we use (4.4a) and (4.7) to deduce that while invoking the assumption (A2) for Ψ ′′ and Ψ ′′′ leads to Then, via Young's inequality Meanwhile, by a similar argument, and so collecting the above estimates for the right-hand side of (4.5) we deduce the existence of two positive constants c 1 and c 2 independent of R such that where we have also used the lower bound (4.6). In the case β > 0 we can directly employ Gronwall's inequality to handle the terms on the right-hand side, whereas in the case β = 0 we proceed as follows: if B > 0, we can choose R > max(2c 1 , 2c 2 B ), and if B = 0 then κ > 0 by (A1) and we employ the generalised Poincaré inequality to handle the term c 2 ψ 2 L 2 on the left-hand side after choosing R sufficiently large. Invoking Gronwall's inequality, keeping in mind that ϕ ∈ L 4 (0, T ; H 2 (Ω)), there exists a constant C independent of β such that (4.9) In view of ξ(0) = 0, we note that from (4.3) the initial data v 0 assigned to v satisfies the elliptic equation Testing with v 0 and using Korn's inequality shows that which explains the absence of initial data on the right-hand side of (4.9). Then, recalling the lower bound (4.6) and employing Korn's inequality we have which also implies the uniqueness of solution since the difference of two solutions to the linear system (4.3) satisfies (4.3) with h 1 = h 2 = h 3 = 0. To complete the proof we return to (4.8) to deduce that while if β > 0, from (4.4c) we also have Lastly after passing to the limit in the Galerkin approximation we obtain limit functions (ξ, η, ψ, v) ∈ Y β lin satisfying (4.3) except for the L 2 (0, T ; H 2 n (Ω)) regularity of ξ. This can be obtained from viewing (4.4b) as the variational formulation of the elliptic problem with a right-hand sidef ∈ L 2 (Q), and with the help of elliptic regularity we then infer that ξ ∈ L 2 (0, T ; H 2 n (Ω)). Hence, we have shown that (ξ, η, ψ, v) ∈ Y β lin and this concludes the proof.

Differentiability of the solution operator
In this section we establish the Fréchet differentiability of the solution operator S between suitable Banach spaces, and that the derivative at w * = (w * 1 , w * 2 , w * 3 ) ∈ U ad in direction h = (h 1 , h 2 , h 3 ) ∈ U is the unique solution (ξ, η, ψ, v) obtained from Theorem 4. This is formulated as follows.
Theorem 5. Under (A1)-(A6), for given w * ∈ U ad with (ϕ, µ, σ, u) = S(w * ) ∈ Y β , the control-to-state operator S is Fréchet differentiable at w * when viewed as a mapping from U to X β , where is the unique solution to (4.3) associated to h.
Proof. We denote and aim to show This is done via establishing for functions with a positive constant C independent of (Φ, λ, θ, z) and h. To this end, we recall from Theorem 1 and Theorem 4 that the new variables (Φ, λ, θ, z) ∈ Y β lin satisfy for all ζ ∈ H 1 (Ω) and η ∈ X(Ω) and for a.e. t ∈ (0, T ), where Take note that h 1 does not appear in (4.11) since the state system (1.1) is linear in w * 1 = σ B . We invoke Taylor's theorem with integral remainder for f ∈ W 2,∞ (R): to deduce that and in light of the regularity assumption (A6), as well as (A2) and (2.4), there exists a positive constant C such that Next, we test (4.11a) with Φ, (4.11b) with λ, (4.11c) with M θ and (4.11d) with Kz for positive constants M, K yet to be determined, and upon adding the resulting equations we arrive at Looking at the terms involving z, we see for a positive constant C independent of K. Next, for terms involving θ, we employ the boundedness of σ and w * 3 , and the Lipschitz continuity of h to obtain for a positive constant C independent of M . Next, for the terms involving k in (4.13) we similarly have , and for the terms involving Ψ ′ , where we used (A2), (2.5) and (4.12). Lastly, we tackle the term involving X h . First, we observe with the assumption g ∈ W 2,∞ (R d×d , R) from (A6) that where the fourth order tensor g ′′ (·) is evaluated at Hence, using the boundedness of f (ϕ), g(·), g ′ (·), g ′′ (·), σ and σ h that are independent of h U , we obtain for a positive constant ε > 0 to be determined later, To close the estimate we require an estimate for Φ 2 H 2 , which can be obtained from (4.11b). Using that (Φ, λ, θ, z) ∈ Y β lin we see that By elliptic regularity, there exists a positive constant C 0 independent of (Φ, λ, θ, z) and h, as well as M and K, such that (4.14) Let α be a positive constant such that Multiplying (4.14) with α and adding to (4.13), then employing the estimates for the right-hand side and choosing ε = α 2 , we obtain where the positive constantsĈ appearing on the left-hand side are independent of M and K. Hence, choosing M and K sufficiently large, with Gronwall's inequality and Korn's inequality, as well as Φ(0) = θ(0) = 0, we have where the last inequality comes from the application of (2.5). This completes the proof as (4.10) has been shown.
Proof. We proceed with formal estimates that can be rigorously derived with a standard Galerkin approximation and let us note that in the following positive constants denoted by the symbol C will be independent of the Galerkin parameter. Then, testing ζ = Kp in (4.16a), φ = −Kq and φ = p in (4.16b), φ = Hr in (4.16c) and η = Zs in (4.16d) for some positive constants K, H and Z yet to be determined, we obtain after summing the resulting equalities (4.17) Firstly, we obtain from (4.16b) and elliptic regularity that Then, a short calculation shows that with positive constants C independent of H and Z. Adding (4.18) to (4.17) and substituting the above yields If B > 0, we choose HB > C, otherwise we use the generalised Poincaré inequality with H sufficiently large so that H ∇r 2 L 2 + Hκ r 2 Then, choosing Z sufficiently large so that Zc 0 > C, and then finally K sufficiently large, we obtain via Gronwall's inequality (applied backwards in time) and Korn's inequality that Then, from (4.16a) we infer p t L 2 (0,T ;H 2 n (Ω) ′ ) ≤ C 1 + Ψ ′′ (ϕ) L ∞ (0,T ;L 2 ) q L 2 (Q) + C p L 2 (Q) + C r L 2 (Q) , and if β > 0, a comparison of terms in (4.16c) gives r t L 2 (0,T ;H 1 (Ω) ′ ) ≤ C r L 2 (0,T ;H 1 ) + C p L 2 (Q) + C q L 2 (Q) .
These estimates are sufficient to pass to the limit and deduce the existence of a solution (p, q, r, s) ∈ Z β to (4.15) in the sense that (4.16) is fulfilled. Moreover, as the adjoint system is linear in (p, q, r, s), the difference of any two solutions satisfy (4.16) where the terms involving n(·, ϕ) and ϕ − ϕ Q are absent in f 1 and f 2 . Consequently, we arrive at an analogue of (4.19) for the difference of two solutions where the right-hand side is zero, which in turn leads to uniqueness of solutions.
Combining these two leads to the simplification and (4.2) is then a consequence of (4.22).

Sparsity of non-negative optimal controls
In the medical context, the control variables w 2 = m(t) and w 3 = s(t) should be nonnegative, and so we modify the admissible control subsets U (2) ad and U where w is a fixed positive constant. If γ 1 > 0, from the optimality condition (4.2), substituting y 2 = w * 2 and y 3 = w * 3 , then using the Hilbert projection theorem allows us to infer that w * 1 is the L 2 (Σ)-orthogonal projection of −κr/γ 1 onto the closed and convex subset U (1) ad of L 2 (Σ), leading to the representation formula w * 1 (x, t) = min w 1 (x, t), max w 1 (x, t), − κ γ 1 r(x, t) for a.e. (x, t) ∈ Σ.
Due to the L 1 -regularisation for w 2 and w 3 in the optimal control problem, we can expect the optimal controls w * 2 and w * 3 to vanish on certain parts of the time interval [0, T ]. This is formulated as follows.