Mathematical analysis of a model-constrained inverse problem for the
reconstruction of early states of prostate cancer growth

Elena Beretta1
e-mail<EMAIL_ADDRESS>
Cecilia Cavaterra2
e-mail<EMAIL_ADDRESS>
Matteo Fornoni3
e-mail<EMAIL_ADDRESS>
Guillermo Lorenzo4
e-mail<EMAIL_ADDRESS>
Elisabetta Rocca5
e-mail<EMAIL_ADDRESS>
1Division of Science, NYU Abu Dhabi

2Department of Mathematics “F. Enriques”, University of Milan & IMATI-C.N.R.

3Department of Mathematics “F. Casorati”, University of Pavia

4Health Research Institute of Santiago de Compostela, Spain & Oden Institute
for Computational Engineering and Sciences, The University of Texas at Austin,
USA

5Department of Mathematics “F. Casorati”, University of Pavia & IMATI-C.N.R.

###### Abstract

The availability of cancer measurements over time enables the personalised
assessment of tumour growth and therapeutic response dynamics. However, many
tumours are treated after diagnosis without collecting longitudinal data, and
cancer monitoring protocols may include infrequent measurements. To facilitate
the estimation of disease dynamics and better guide ensuing clinical
decisions, we investigate an inverse problem enabling the reconstruction of
earlier tumour states by using a single spatial tumour dataset and a
biomathematical model describing disease dynamics. We focus on prostate
cancer, since aggressive cases of this disease are usually treated after
diagnosis. We describe tumour dynamics with a phase-field model driven by a
generic nutrient ruled by reaction-diffusion dynamics. The model is completed
with another reaction-diffusion equation for the local production of prostate-
specific antigen, which is a key prostate cancer biomarker. We first improve
previous well-posedness results by further showing that the solution operator
is continuously Fréchet differentiable. We then analyse the backward inverse
problem concerning the reconstruction of earlier tumour states starting from
measurements of the model variables at the final time. Since this problem is
severely ill-posed, only very weak conditional stability of logarithmic type
can be recovered from the terminal data. However, by restricting the unknowns
to a compact subset of a finite-dimensional subspace, we can derive an optimal
Lipschitz stability estimate.

Key words: phase field; non-linear parabolic system; inverse problems; well-
posedness; mathematical oncology; prostate cancer.

AMS (MOS) Subject Classification: 35K51, 35K58 35R30, 35Q92, 92C50

##  1 Introduction

The mathematical modelling of cancer growth and therapeutic response has
enabled the investigation of the biophysical mechanisms underlying the
development and treatment of these diseases, as well as the prediction of
personalised outcomes to guide clinical decision-making [55, 41]. Personalised
tumour forecasting requires patient-specific measurements of cancer growth and
treatment response, such as biomarker, imaging, histopathological and omics
data [41, 37]. These measurements enable the calibration of model parameters,
the definition of the tumour and host tissue geometry, and the initialisation
of model variables (e.g., those describing tumour and other cell species as
well as key substances controlling their dynamics) [55, 41, 37]. Towards these
ends, data are usually required at two or more time points during the
monitoring of the disease before, during and after treatment [41, 40, 29, 52].
The collection of these data over time contributes to early characterise
crucial clinical endpoints (e.g., progression to higher-risk disease,
treatment failure, survival) by both analysing the relative change of these
measurements [25, 49] as well as by using them to obtain more accurate tumour
forecasts [41, 40, 29, 52]. However, depending on the natural history and
specific clinical management protocols for each type of cancer, it may not be
possible to obtain tumour measurements until the tumour has developed
sufficiently to produce symptoms or be detected with standard-of-care
screening methods [41, 37, 12, 40]. Furthermore, tumour monitoring strategies
may not enable to assess the tumour status as often as necessary to capture
tumour growth and treatment response in detail [41, 37, 12, 40]. These
situations may lead to limited accuracy in the estimation of disease prognosis
and result in treatment excesses and deficiencies that can respectively affect
the patients’ quality of life and life expectancy [26, 56, 46].

Figure 1: Reconstruction of earlier prostate cancer states for clinical
decision-making. This figure illustrates the potential application of
reconstructed prostate cancer dynamics from a diagnostic medical imaging
dataset (e.g., MRI), while using a mathematical model describing prostate
cancer growth. The first step (1) consists of estimating an earlier tumour
state for the tumour within a clinically reasonable time, for example, six to
twelve months in the context of monitoring prostate cancer or decision of
first-line treatment. This is the step in which we focus on the present work.
Then, a second step (2) leverages the reconstructed tumour growth dynamics to
make a forecast within comparable time horizons to the one used in the
reconstruction of earlier states. In the third and final step (3), these
predictions are analysed to inform clinical decisions, for instance, whether
to schedule another medical imaging scan to monitor lower-risk disease or to
prescribe immediate treatment due to anticipated worsening of the disease in
the short-term future (e.g., surgery, radiotherapy, or drug-based therapies).

After the acquisition of a medical imaging measurement to characterise the
tumour morphology (e.g., at diagnosis, after treatment), an inverse problem
leveraging an adequate mathematical model to represent tumour dynamics and
aiming at identifying the tumour status at an earlier date can shed light on
the patient-specific tumour dynamics [41, 48, 12]. This reconstruction can
then be used to better estimate potential clinical outcomes that guide crucial
decision-making by leveraging tumour forecasts or model-based biomarkers.
Figure 1 illustrates the overall procedure whereby a spatial measurement of a
tumour obtained _via_ medical imaging at the time of diagnosis is used to
reconstruct an earlier tumour state and better inform the ensuing management
of the disease. Indeed, the determination of the spatial region where a tumour
originates can be of clinical relevance for some tumour types, such as brain
and prostate cancers [4, 35]. Additionally, the need to estimate earlier
states of tumour development can also arise in a pre-clinical context. For
example, it may be necessary to determine the configuration of an _in vitro_
cell colony or the tumour morphology in an _in vivo_ animal model before the
onset of experimental measurements or due to sparse data collection with
important tumour dynamics between consecutive data points.

In this work, we analyse a model-constrained inverse problem aiming to
reconstruct earlier tumour configurations from a spatial measurement
collected, for instance, _via_ medical imaging at the time of diagnosis. We
focus on newly detected prostate cancer. The diagnosis of these tumours relies
on a multiparametric magnetic resonance imaging (MRI) scan, which is often
motivated by increasing values of the blood biomarker prostate-specific
antigen (PSA; the most common biomarker in prostate cancer management) and
guides the ensuing biopsy procedure to confirm the disease histopathologically
[45, 40]. The first-line clinical management options for prostate cancer
usually include active surveillance for indolent, lower-risk cases and radical
treatment (e.g., surgery, radiotherapy with/without androgen deprivation
therapy) for aggressive, higher-risk patients [45]. Thus, the accurate
identification of the clinical risk of prostate cancer is a critical need of
utmost importance to select an adequate clinical management strategy. Beyond
the diagnostic PSA, MRI, and biopsy data, estimating the recent dynamics of
prostate cancer growth before diagnosis can further inform the treating
physicians about the best course of management for each particular patient
[25, 40]. Towards this end, we employ a mathematical model of prostate cancer
growth based on the phase-field method, which has been studied both from an
analytical and a numerical point of view in [13] and [16]. This PDE modelling
approach is a common spatiotemporal continuous formulation to describe the
time evolution of the location and geometry of a tumour, as well as to address
related optimal control problems [11, 14, 15, 20, 21, 22, 24, 23, 42, 44, 43,
51, 53]. Hence, we define a continuous phase field $\varphi$ that takes values
$\varphi\approx 0$ in the host tissue, $\varphi\approx 1$ in the tumour, and
intermediate values in a small neighbourhood of the tumour-healthy tissue
interface, which exhibits a smooth yet steep profile. The tumour phase field
dynamics depends on the local availability of a critical nutrient (e.g.,
oxygen or glucose), whose concentration is denoted by $\sigma$ and obeys a
reaction-diffusion equation. Additionally, the model includes another
reaction-diffusion equation to represent the local production of tissue PSA
$p$, which represents the PSA leaked to the bloodstream per unit volume of
prostatic tissue [42]. The complete prostate cancer model considered herein
reads as follows:

$\displaystyle\partial_{t}\varphi-\lambda\Delta\varphi+F^{\prime}(\varphi)-m(\sigma)\mathbbm{h}^{\prime}(\varphi)=0$
in $Q_{T}$, (1.1)
$\displaystyle\partial_{t}\sigma-\eta\Delta\sigma=S_{h}(1-\varphi)+S_{c}\varphi-\left(\gamma_{h}(1-\varphi)+\gamma_{c}\varphi\right)\sigma\qquad$
in $Q_{T}$, (1.2) $\displaystyle\partial_{t}p-D\Delta
p+\gamma_{p}p=\alpha_{h}(1-\varphi)+\alpha_{c}\varphi$ in $Q_{T}$, (1.3)
$\displaystyle\varphi=0,\quad\partial_{\boldsymbol{n}}\sigma=0,\quad\partial_{\boldsymbol{n}}p=0$
on $\Sigma_{T}$, (1.4)
$\displaystyle\varphi(0)=\varphi_{0},\quad\sigma(0)=\sigma_{0},\quad
p(0)=p_{0}$ in $\Omega$, (1.5)

where $\Omega\subset\mathbb{R}^{N}$, $N=2,3$, is a sufficiently regular open
and bounded domain with outward unit normal vector $\boldsymbol{n}$, $T>0$ is
the final time (or time horizon), $Q_{T}:=\Omega\times(0,T)$, and
$\Sigma_{T}=\partial\Omega\times(0,T)$. For the inverse problem addressed
herein, $t=T$ represents the time of the diagnostic imaging scan while $t=0$
is a previous timepoint at which we want to reconstruct the tumour phase
field.

The diffusion coefficient of the tumour phase in Eq. (1.1) can be defined as
$\lambda=M\ell^{2}$, where both $M$ and $\ell$ are positive real constants
denoting the tumour mobility and interface width [53]. The non-linear
functions $F$ and $\mathbbm{h}$ are explicitly defined as:

$\displaystyle
F(\varphi)=M\varphi^{2}(1-\varphi)^{2},\qquad\mathbbm{h}(\varphi)=M\varphi^{2}(3-2\varphi),\qquad\text{with
}M>0.$

The non-convex function $F(\varphi)$ is a double-well potential, which is
typical in phase field modelling and allows for the coexistence of the tumour
($\varphi\approx 1$) and healthy ($\varphi\approx 0$) tissues. The
interpolation function $\mathbbm{h}(\varphi)$ is also common in non-conserved
dynamics of phase fields and it verifies the properties $\mathbbm{h}(0)=0$ and
$\mathbbm{h}(1)=1$ as well as
$\mathbbm{h}^{\prime}(0)=\mathbbm{h}^{\prime}(1)=0$, where
$\mathbbm{h}^{\prime}$ denotes the derivative of $\mathbbm{h}$. The function
$m(\sigma)$ is usually called tilting function and, in our model, is defined
as

$m(\sigma)=m_{ref}\left(\frac{\rho+A}{2}+\frac{\rho-A}{\pi}\arctan\left(\frac{\sigma-\sigma_{l}}{\sigma_{r}}\right)\right),$
(1.6)

where $m_{ref}$ is a positive constant that scales the strength of the tilting
function within our phase field framework, while $\rho$ and $A$ are constant
proliferation and apoptosis indices, respectively. These non-dimensional
parameters are defined upon the proliferation and apoptosis rates in tumour
tissue as follows:

$\rho=\frac{K_{\rho}}{\overline{K_{\rho}}},\quad\text{and}\quad
A=-\frac{K_{A}}{\overline{K_{A}}},$

where $K_{\rho}$ is the proliferation rate of tumour cells, $K_{A}$ is the
apoptosis rate of tumour cells, $\overline{K_{\rho}}$ is a scaling reference
value for the proliferation rate, and $\overline{K_{A}}$ is a scaling
reference value for the apoptosis rate. Thus, according to the definition of
$m(\sigma)$ provided above, this function represents the nutrient-dependent
net proliferation rate of the tumour. Additionally, the constants $\sigma_{r}$
and $\sigma_{l}$ in Eq. (1.6) denote a reference and a threshold value for the
nutrient concentration[53]. Therefore, nutrient concentrations lower than
$\sigma_{l}$ render the healthy tissue energetically more favorable than the
tumour tissue and vice-versa. Moreover, when $|m(\sigma)|<1/3$, the function
$G(\varphi,\sigma)=F(\varphi)-m(\sigma)\mathbbm{h}(\varphi)$ is still a
double-well potential presenting two local minima at $\varphi=0$ and
$\varphi=1$ [13]. Within this range, low nutrient concentrations result in a
lower energy level (i.e., the value of $G$) in the healthy tissue
($\varphi=0$) than in the tumoral tissue ($\varphi=1$). The opposite is true
if we consider high values of the nutrient concentration $\sigma$. The
described behaviours of $m(\sigma)$ and $G(\varphi,\sigma)$ are further
illustrated in Figure 2.

(a)

(b)

(c)

Figure 2: Nonlinear terms in the phase field model of prostate cancer growth.
Panel (a) depicts the tilting function $m(\sigma)$, which represents the net
proliferation, with respect to the local nutrient $\sigma$. This plot shows
that proliferation is supported in regions with higher nutrient
concentrations, while tumour cell death dominates in areas where the nutrient
is scarcer. Panel (b) shows the potential
$G\left(\varphi,\sigma\right)=F(\varphi)-m(\sigma)\mathbbm{h}(\varphi)$ with
respect to the values of the tumour phase field $\varphi$ in three situations:
a high nutrient availability that fosters proliferation ($m(\sigma)>0)$,
limited nutrient availability that inhibits proliferation ($m(\sigma)=0$), and
a scarce nutrient concentration that promotes tumour cell death
($m(\sigma)<0$). For the same three scenarios, panel (c) further provides the
derivative $\partial
G\left(\varphi,\sigma\right)/\partial\varphi=F^{\prime}(\varphi)-m(\sigma)\mathbbm{h}^{\prime}(\varphi)$,
which explicitly drives the spatiotemporal dynamics of the tumour phase field
in Eq. (2.1). The plots in panels (b) and (c) show that a high nutrient
concentration energetically favours tumour growth within our modelling
framework, whereas a low nutrient environment energetically impedes it.

In the nutrient dynamics equation (Eq. (1.2)), $\eta$ denotes the diffusion
coefficient of the nutrient, $S_{h}$ and $S_{c}$ stand for the non-negative
nutrient supply rates in the healthy and tumour tissue, and $\gamma_{h}$ and
$\gamma_{c}$ are positive constants that represent the nutrient uptake rate in
the healthy and cancerous tissue. Additionally, the parameters governing the
tissue PSA equation (Eq. 1.3) are the diffusion coefficient $D$, the rates of
production of tissue PSA by the healthy tissue $\alpha_{h}$ and the tumour
$\alpha_{c}$, and the natural decay rate $\gamma_{p}$. To conclude the
exposition of the model, we choose standard no-flux boundary conditions for
$\sigma$ and $p$, while we set zero-valued Dirichlet boundary conditions for
$\varphi$. This last assumption implies that the tumour is confined inside the
prostate (recall indeed that $\varphi\approx 0$ represents the healthy phase),
which is the case of the majority of newly-diagnosed prostate cancer cases
[45]. Furthermore, the only difference with respect to the original model
presented in [13] is that we have neglected the terms related to cytotoxic and
antiangiogenic treatments. The rationale for this choice is that here we focus
on reconstructing early stages of prostate cancer growth before the time of
first diagnosis using imaging data and, thus, before the administration of any
treatment [40, 45]. The interested reader is referred to [13, 16] for a more
detailed presentation of the model equations, parameters, and alternative
definitions (e.g., further mechanisms to include in the model, spatially-
varying parameters, nonlinearities).

Mathematically, the problem we want to address is a backward inverse problem
for a nonlinear system of parabolic equations. Such kind of problems arise any
time a diffusion process has to be reversed thereby governing a multitude of
applications [5, 34, 48, 18, 2]. To be precise, we call
$\mathcal{R}:(\varphi_{0},\sigma_{0},p_{0})\mapsto(\varphi(T),\sigma(T),p(T))$
the nonlinear solution operator which associates to any initial data
$(\varphi_{0},\sigma_{0},p_{0})$ the value of the solution at the final time
$(\varphi(T),\sigma(T),p(T))$. Then, given some measurements
$\varphi_{\text{meas}},\sigma_{\text{meas}},p_{\text{meas}}:\Omega\to\mathbb{R}$,
we can formulate our inverse problem as that of finding initial data
$\varphi_{0}$, $\sigma_{0}$, $p_{0}$ such that the solution
$(\varphi,\sigma,p)$ to (1.1)–(1.5) satisfies
$\varphi(T)=\varphi_{\text{meas}}$, $\sigma(T)=\sigma_{\text{meas}}$ and
$p(T)=p_{\text{meas}}$. The backward inverse problem, even in the linear case,
is well-known to be severely ill-posed with very weak conditional logarithmic
stability from the data [47, 39, 33]. Therefore, as a preliminary step, it is
of paramount importance to regularise the problem, by establishing physically
relevant a priori assumptions on the unknown initial data that lead to better,
possibly Lipschitz, dependence of the initial data on the measurements. In
turn, such stability estimates give local convergence of iterative
reconstruction methods [17]. This strategy has been successfully employed in
addressing other non-linear severely ill-posed problems, for example in the
context of the conductivity inverse problem, see for example [3, 9, 6]. The
key idea is that, by restricting the unknowns to a compact subset of a finite-
dimensional subspace, we can derive Lipschitz dependence, via a quantitative
version of the inverse map theorem in Banach spaces [7, Proposition 5]. To
accomplish this, we need to prove the injectivity and continuity of the
inverse map, as well as Fréchet differentiability with continuity of the
forward map. Indeed, even though the system is non-linear, we can adapt
logarithmic convexity methods to prove injectivity and a conditional stability
estimate for the inverse map for any positive time $t>0$ (Proposition 3.1).
Additionally, through a more refined analysis, we are also able to recover a
logarithmic stability estimate for the reconstruction of the initial data
(Theorem 3.5). We just mention that the log-convexity method was already
employed for example in [28] to derive stability estimates for linear
parabolic partial differential equations and in [38] in the case of non-linear
parabolic equations. Carleman estimates for parabolic equations can be also
employed to derive uniqueness and stability estimates for the backward inverse
problem, see for instance [54, 31, 32, 50]. Then, to prove the additional
regularity of $\mathcal{R}$ required by [7, Proposition 5], we show finer
continuous dependence estimates on the solutions of the system (1.1)–(1.5),
that allow us to further get that $\mathcal{R}$ is $\mathcal{C}^{1}$ (Theorem
2.5). As a consequence, we can derive the sought quantitative Lipschitz
stability estimate (Theorem 3.10), under the a priori assumption that the
initial data lie in a compact subset of a finite-dimensional subspace. In
particular, our explicit computation of the stability constants confirms that,
as expected, the inverse problem is severely ill-posed. Indeed, such constants
blow up exponentially with the dimension of the finite-dimensional subspace
and doubly exponentially with the final time $T$. However, such a Lipschitz
stability estimate paves the way for the use of iterative locally convergent
reconstruction algorithms, such as the Landweber iteration method [17, 36].
Due to the challenging nature of the ill-posed inverse problem, we plan to
address an in-depth numerical study in a forthcoming paper, aiming to explore
the possibility of an efficient reconstruction to be eventually used in
clinical procedures.

The paper is organised as follows. In Section 2 we recall the well-posedness
results shown in [13] and then prove that the solution operator is
continuously Fréchet differentiable. Section 3 is devoted to the analysis of
the inverse problem of reconstructing the initial data. Here we prove our main
stability results in Theorem 3.5 and Theorem 3.10. Section 4 finally
summarises the obtained results and explores future directions for the
numerical study of the problem.

##  2 Well-posedness and differentiability of the solution operator

We first introduce some notation that will be used throughout the paper. Let
$\Omega\subset\mathbb{R}^{N}$, $N=2,3$, be an open and bounded domain with
$\mathcal{C}^{2}$ boundary and outward unit normal vector $\boldsymbol{n}$.
Let $T>0$ be the final time, and denote $Q_{t}:=\Omega\times(0,t)$ and
$\Sigma_{t}=\partial\Omega\times(0,t)$ for any $t\in(0,T]$. We define the
spaces:

$H=L^{2}(\Omega),\quad V_{0}=H^{1}_{0}(\Omega),\quad
V_{0}^{*}=(H^{1}_{0}(\Omega))^{*},\quad V=H^{1}(\Omega),\quad
V^{*}=(H^{1}(\Omega))^{*},$ $W_{0}=H^{2}(\Omega)\cap H^{1}_{0}(\Omega),\quad
W=\left\\{u\in H^{2}(\Omega)\mid\partial_{\boldsymbol{n}}u=0\right\\}.$

By standard results, we know that, if $H$ is identified with its dual, the
following compact and dense embeddings hold:

$W_{0}\hookrightarrow V_{0}\hookrightarrow H\hookrightarrow
V_{0}^{*}\quad\text{and}\quad W\hookrightarrow V\hookrightarrow
H\hookrightarrow V^{*}.$

Recall that, by elliptic regularity, we can use the equivalent norms:

$\lVert u\rVert^{2}_{W_{0}}=\lVert u\rVert^{2}_{H}+\lVert\Delta
u\rVert^{2}_{H},\quad\lVert u\rVert^{2}_{W}=\lVert
u\rVert^{2}_{H}+\lVert\Delta u\rVert^{2}_{H}.$

In some cases, we will denote for simplicity:

$\mathbb{H}=H\times H\times H,\quad\mathbb{V}=V_{0}\times V\times V.$

By using standard Bochner theory for spaces of time-dependent functions with
values on Banach spaces, we also introduce the spaces in which weak and strong
solutions of system (2.1)–(2.5) will be respectively defined:

$\displaystyle\mathbb{Y}_{0}$
$\displaystyle=H^{1}(0,T;V_{0}^{*})\cap\mathcal{C}^{0}([0,T];H)\cap
L^{2}(0,T;V_{0}),$ $\displaystyle\mathbb{Y}$
$\displaystyle=H^{1}(0,T;V^{*})\cap\mathcal{C}^{0}([0,T];H)\cap L^{2}(0,T;V),$
$\displaystyle\mathbb{X}_{0}$
$\displaystyle=H^{1}(0,T;H)\cap\mathcal{C}^{0}([0,T];V_{0})\cap
L^{2}(0,T;W_{0}),$ $\displaystyle\mathbb{X}$
$\displaystyle=H^{1}(0,T;H)\cap\mathcal{C}^{0}([0,T];V)\cap L^{2}(0,T;W).$

Regarding the parameters of the system, we assume the following hypotheses:

  1. A1.

$\lambda,\eta,\gamma_{h},\gamma_{c},S_{h},S_{c},D,\gamma_{p},\alpha_{h},\alpha_{c}>0$.

  2. A2.

$F(s)=Ms^{2}(1-s)^{2}$ and $\mathbbm{h}(s)=Ms^{2}(3-2s)$, with $M>0$, for any
$s\in\mathbb{R}$. In particular, we observe that
$F,\mathbbm{h}\in\mathcal{C}^{3}(\mathbb{R})$.

  3. A3.

$m(s)=m_{\text{ref}}\left(\frac{\rho+A}{2}+\frac{\rho-A}{2}\arctan\left(\frac{s-\sigma_{l}}{\sigma_{r}}\right)\right)$,
with $m_{\text{ref}},\rho,A,\sigma_{l},\sigma_{r}>0$. In particular, we
observe that $m$ and $m^{\prime}$ are Lipschitz continuous on $\mathbb{R}$ and
$m,m^{\prime},m^{\prime\prime}\in L^{\infty}(\mathbb{R})$.

In what follows, we will extensively use the symbol $C>0$, which may also
change from line to line, to denote positive constants depending only on the
fixed parameters of the system and possibly on $\Omega$ and $T$. In some
cases, we will use a subscript to highlight some particular dependence of
these constants. Finally, to simplify a bit the analysis of the source terms
in the model, we rewrite it in the following way:

$\displaystyle\partial_{t}\varphi-\lambda\Delta\varphi+F^{\prime}(\varphi)-m(\sigma)\mathbbm{h}^{\prime}(\varphi)=0$
in $Q_{T}$, (2.1)
$\displaystyle\partial_{t}\sigma-\eta\Delta\sigma=S_{h}+S_{ch}\varphi-\gamma_{h}\sigma-\gamma_{ch}\sigma\varphi\qquad$
in $Q_{T}$, (2.2) $\displaystyle\partial_{t}p-D\Delta
p+\gamma_{p}p=\alpha_{h}+\alpha_{ch}\varphi$ in $Q_{T}$, (2.3)
$\displaystyle\varphi=0,\quad\partial_{\boldsymbol{n}}\sigma=\partial_{\boldsymbol{n}}p=0$
in $\Sigma_{T}$, (2.4)
$\displaystyle\varphi(0)=\varphi_{0},\quad\sigma(0)=\sigma_{0},\quad
p(0)=p_{0}$ in $\Omega$, (2.5)

where we call $\gamma_{ch}:=\gamma_{c}-\gamma_{h}$, $S_{ch}:=S_{c}-S_{h}$ and
$\alpha_{ch}:=\alpha_{c}-\alpha_{h}$.

Our first step consists of the analysis of the forward system. We start by
recalling the following result about the well-posedness of system (2.1)–(2.5),
proved in [13, Theorem 3.2].

######  Theorem 2.1.

Assume hypotheses A1–A3 and let

$(\varphi_{0},\sigma_{0},p_{0})\in\mathbb{H},\quad\text{and}\quad
0\leq\varphi_{0}(x)\leq 1\quad\text{for a.e. }x\in\Omega.$ (2.6)

Then, there exists a unique weak solution

$(\varphi,\sigma,p)\in\mathbb{Y}_{0}\times\mathbb{Y}\times\mathbb{Y},$

with

$0\leq\varphi(x,t)\leq 1\quad\text{a.e. }(x,t)\in Q_{T},$ (2.7)

which solves system (2.1)–(2.5) in variational formulation, namely

$\displaystyle\langle\partial_{t}\varphi,v\rangle_{V_{0}}+\lambda(\nabla\varphi,\nabla
v)_{H}+(F^{\prime}(\varphi)-m(\sigma)\mathbbm{h}^{\prime}(\varphi),v)_{H}=0$
$\displaystyle\hbox{a.e. in $(0,T)$, $\forall v\in V_{0}$},$
$\displaystyle\langle\partial_{t}\sigma,v\rangle_{V}+\eta(\nabla\sigma,\nabla
v)_{H}+(\gamma_{h}\sigma+\gamma_{ch}\sigma\varphi-
S_{h}-S_{ch}\varphi,v)_{H}=0$ $\displaystyle\hbox{a.e. in $(0,T)$, $\forall
v\in V$},$ $\displaystyle\langle\partial_{t}p,v\rangle_{V}+D(\nabla p,\nabla
v)_{H}+(\gamma_{p}p-\alpha_{h}-\alpha_{ch}\varphi,v)_{H}=0$
$\displaystyle\hbox{a.e. in $(0,T)$, $\forall v\in V$},$

with $\varphi(0)=\varphi_{0}$, $\sigma(0)=\sigma_{0}$ and $p(0)=p_{0}$, and
satisfies

$\lVert\varphi\rVert^{2}_{\mathbb{Y}_{0}}+\lVert\sigma\rVert^{2}_{\mathbb{Y}}+\lVert
p\rVert^{2}_{\mathbb{Y}}\leq
C\left(1+\lVert\varphi_{0}\rVert^{2}_{H}+\lVert\sigma_{0}\rVert^{2}_{H}+\lVert
p_{0}\rVert^{2}_{H}\right),$ (2.8)

for some $C>0$ depending only on the parameters of the system. Moreover, if

$\begin{split}\sigma_{0},\,p_{0}\in
L^{\infty}(\Omega),\quad\text{and}\quad\sigma_{0}(x)\geq 0,\,p_{0}(x)\geq
0\text{ for a.e. }x\in\Omega,\end{split}$ (2.9)

then also

$\sigma,p\in L^{\infty}(Q_{T}),\quad\text{and}\quad\sigma(x,t)\geq
0,\,p(x,t)\geq 0\text{ for a.e. }(x,t)\in Q_{T},$

with

$\lVert\sigma\rVert_{L^{\infty}(Q_{T})}\leq
C(\lVert\sigma_{0}\rVert_{L^{\infty}(\Omega)}+1),\quad\lVert
p\rVert_{L^{\infty}(Q_{T})}\leq C(\lVert p_{0}\rVert_{L^{\infty}(\Omega)}+1).$

Furthermore, the solution is continuous with respect to the initial data, that
is, taking two solutions $(\varphi_{i},\sigma_{i},p_{i})$ corresponding to
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})$, $i=1,2$, we have

$\begin{split}&\lVert(\varphi_{1}-\varphi_{2})(t)\rVert^{2}_{H}+\lVert(\sigma_{1}-\sigma_{2})(t)\rVert^{2}_{H}+\lVert(p_{1}-p_{2})(t)\rVert^{2}_{H}\\\
&\quad+\lVert\varphi_{1}-\varphi_{2}\rVert^{2}_{L^{2}(0,T;V_{0})}+\lVert\sigma_{1}-\sigma_{2}\rVert^{2}_{L^{2}(0,T;V)}+\lVert
p_{1}-p_{2}\rVert^{2}_{L^{2}(0,T;V)}\\\ &\qquad\leq
C\left(\lVert\varphi_{0}^{1}-\varphi_{0}^{2}\rVert^{2}_{H}+\lVert\sigma_{0}^{1}-\sigma_{0}^{2}\rVert^{2}_{H}+\lVert
p_{0}^{1}-p_{0}^{2}\rVert^{2}_{H}\right),\end{split}$ (2.10)

for any $t\in[0,T]$. Finally, if
$(\varphi_{0},\sigma_{0},p_{0})\in\mathbb{V}$, the solution enjoys the higher
regularity
$(\varphi,\sigma,p)\in\mathbb{X}_{0}\times\mathbb{X}\times\mathbb{X}$, solves
the system (2.1)–(2.5) almost everywhere and satisfies the estimate

$\lVert\varphi\rVert^{2}_{\mathbb{X}_{0}}+\lVert\sigma\rVert^{2}_{\mathbb{X}}+\lVert
p\rVert^{2}_{\mathbb{X}}\leq
C\left(1+\lVert\varphi_{0}\rVert^{2}_{V_{0}}+\lVert\sigma_{0}\rVert^{2}_{V}+\lVert
p_{0}\rVert^{2}_{V}\right).$ (2.11)

######  Remark 2.2.

Due to our model, hypotheses A2 and A3 are very specific on the non-
linearities that we consider. However, from an analytical point of view, one
could consider more general expressions of $F$, $\mathbbm{h}$ and $m$, as long
as the results of Theorem 2.1 still hold. Hence, the forthcoming analysis can
easily be extended to a more general framework.

We additionally prove a stronger continuous dependence estimate, which will be
useful when proving the Fréchet-differentiability of the solution operator.

######  Proposition 2.3.

Assume hypotheses A1–A3 and let
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})\in\mathbb{V}$, $i=1,2$, be two
sets of initial data, satisfying also (2.6) and (2.9). Let
$(\varphi_{i},\sigma_{i},p_{i})$, $i=1,2$, be the corresponding strong
solutions of (2.1)–(2.5). Then, there exists a constant $C>0$, depending only
on the parameters of the system and on
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})$, $i=1,2$, but not on their
difference, such that

$\begin{split}&\lVert(\varphi_{1}-\varphi_{2})(t)\rVert^{2}_{V_{0}}+\lVert(\sigma_{1}-\sigma_{2})(t)\rVert^{2}_{V}+\lVert(p_{1}-p_{2})(t)\rVert^{2}_{V}\\\
&\quad+\lVert\varphi_{1}-\varphi_{2}\rVert^{2}_{L^{2}(0,T;W_{0})}+\lVert\sigma_{1}-\sigma_{2}\rVert^{2}_{L^{2}(0,T;W)}+\lVert
p_{1}-p_{2}\rVert^{2}_{L^{2}(0,T;W)}\\\ &\qquad\leq
C\left(\lVert\varphi_{0}^{1}-\varphi_{0}^{2}\rVert^{2}_{V_{0}}+\lVert\sigma_{0}^{1}-\sigma_{0}^{2}\rVert^{2}_{V}+\lVert
p_{0}^{1}-p_{0}^{2}\rVert^{2}_{V}\right),\quad\forall\,t\in[0,T].\end{split}$
(2.12)

###### Proof.

Let us define the differences

$\varphi:=\varphi_{1}-\varphi_{2},\quad\sigma:=\sigma_{1}-\sigma_{2},\quad
p:=p_{1}-p_{2}.$

Then, up to adding and subtracting some terms, they satisfy the following
system of equations:

$\displaystyle\partial_{t}\varphi-\lambda\Delta\varphi+(F^{\prime}(\varphi_{1})-F^{\prime}(\varphi_{2}))$
$\displaystyle\qquad-m(\sigma_{1})(\mathbbm{h}^{\prime}(\varphi_{1})-\mathbbm{h}^{\prime}(\varphi_{2}))-(m(\sigma_{1})-m(\sigma_{2}))\mathbbm{h}^{\prime}(\varphi_{2})=0$
in $Q_{T}$, (2.13)
$\displaystyle\partial_{t}\sigma-\eta\Delta\sigma+\gamma_{h}\sigma+\gamma_{ch}\sigma_{1}(\varphi_{1}-\varphi_{2})+\gamma_{ch}(\sigma_{1}-\sigma_{2})\varphi_{2}=S_{ch}\varphi\qquad$
in $Q_{T}$, (2.14) $\displaystyle\partial_{t}p-D\Delta
p+\gamma_{p}p=\alpha_{ch}\varphi$ in $Q_{T}$, (2.15)
$\displaystyle\varphi=0,\quad\partial_{\boldsymbol{n}}\sigma=\partial_{\boldsymbol{n}}p=0$
in $\Sigma_{T}$, (2.16)
$\displaystyle\varphi(0)=\varphi_{0}^{1}-\varphi_{0}^{2},\quad\sigma(0)=\sigma_{0}^{1}-\sigma_{0}^{2},\quad
p(0)=p_{0}^{1}-p_{0}^{2}$ in $\Omega$. (2.17)

Note that, being the initial data in $\mathbb{V}$, by Theorem 2.1 the
equalities above have to be intended in a strong sense. Therefore, we are
allowed to test (2.13) by $-\Delta\varphi$, (2.14) by $-\Delta\sigma$, (2.15)
by $-\Delta p$ and sum them up, to obtain:

$\begin{split}&\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\lVert\nabla\varphi\rVert^{2}_{H}+\lVert\nabla\sigma\rVert^{2}_{H}+\lVert\nabla
p\rVert^{2}_{H}\right)+\lambda\lVert\Delta\varphi\rVert^{2}_{H}+\eta\lVert\Delta\sigma\rVert^{2}_{H}+D\lVert\Delta
p\rVert^{2}_{H}\\\
&\quad=-\gamma_{h}\lVert\nabla\sigma\rVert^{2}_{H}-\gamma_{p}\lVert\nabla
p\rVert^{2}_{H}+(F^{\prime}(\varphi_{1})-F^{\prime}(\varphi_{2}),\Delta\varphi)_{H}\\\
&\qquad-(m(\sigma_{1})(\mathbbm{h}^{\prime}(\varphi_{1})-\mathbbm{h}^{\prime}(\varphi_{2})),\Delta\varphi)_{H}-((m(\sigma_{1})-m(\sigma_{2}))\mathbbm{h}^{\prime}(\varphi_{2}),\Delta\varphi)_{H}\\\
&\qquad+\gamma_{ch}(\sigma_{1}(\varphi_{1}-\varphi_{2}),\Delta\sigma)_{H}+\gamma_{ch}((\sigma_{1}-\sigma_{2})\varphi_{2},\Delta\sigma)_{H}\\\
&\qquad-S_{ch}(\varphi,\Delta\sigma)_{H}-\alpha_{ch}(\varphi,\Delta
p)_{H}.\end{split}$ (2.18)

Then, we can easily estimate the terms on the right-hand side by using Cauchy-
Schwarz and Young’s inequalities, together with the local Lipschitz continuity
of $F^{\prime}$, $\mathbbm{h}^{\prime}$, $m$ and the fact that
$\varphi_{i},\sigma_{i}\in L^{\infty}(Q_{T})$ for $i=1,2$. Indeed, we have
that:

$\displaystyle(F^{\prime}(\varphi_{1})-F^{\prime}(\varphi_{2}),\Delta\varphi)_{H}-(m(\sigma_{1})(\mathbbm{h}^{\prime}(\varphi_{1})-\mathbbm{h}^{\prime}(\varphi_{2})),\Delta\varphi)_{H}-((m(\sigma_{1})-m(\sigma_{2}))\mathbbm{h}^{\prime}(\varphi_{2}),\Delta\varphi)_{H}$
$\displaystyle\qquad+\gamma_{ch}(\sigma_{1}(\varphi_{1}-\varphi_{2}),\Delta\sigma)_{H}+\gamma_{ch}((\sigma_{1}-\sigma_{2})\varphi_{2},\Delta\sigma)_{H}-S_{ch}(\varphi,\Delta\sigma)_{H}-\alpha_{ch}(\varphi,\Delta
p)_{H}$
$\displaystyle\quad\leq\frac{\lambda}{2}\lVert\Delta\varphi\rVert^{2}_{H}+\frac{\eta}{2}\lVert\Delta\sigma\rVert^{2}_{H}+\frac{D}{2}\lVert\Delta
p\rVert^{2}_{H}+C\lVert\varphi\rVert^{2}_{H}+C\lVert\sigma\rVert^{2}_{H}.$

Hence, by also integrating on $(0,t)$, for any $t\in(0,T)$, we get that

$\displaystyle\frac{1}{2}\left(\lVert\nabla\varphi(t)\rVert^{2}_{H}+\lVert\nabla\sigma(t)\rVert^{2}_{H}+\lVert\nabla
p(t)\rVert^{2}_{H}\right)+\frac{\lambda}{2}\int_{0}^{t}\lVert\Delta\varphi\rVert^{2}_{H}\,\mathrm{d}s+\frac{\eta}{2}\int_{0}^{t}\lVert\Delta\sigma\rVert^{2}_{H}\,\mathrm{d}s+\frac{D}{2}\int_{0}^{t}\lVert\Delta
p\rVert^{2}_{H}\,\mathrm{d}s$
$\displaystyle\quad\leq\frac{1}{2}\left(\lVert\varphi_{0}^{1}-\varphi_{0}^{2}\rVert^{2}_{V_{0}}+\lVert\sigma_{0}^{1}-\sigma_{0}^{2}\rVert^{2}_{V}+\lVert
p_{0}^{1}-p_{0}^{2}\rVert^{2}_{V}\right)+C\int_{0}^{T}\lVert\varphi\rVert^{2}_{V_{0}}\,\mathrm{d}s+C\int_{0}^{T}\lVert\sigma\rVert^{2}_{V}\,\mathrm{d}s,$

which by Gronwall’s inequality and (2.10) implies (2.12). Proposition 2.3 is
proved. ∎

Given the results of Theorem 2.1 and Proposition 2.3, we introduce the set of
admissible initial data as follows:

$\mathcal{I}_{\text{ad}}=\left\\{(\varphi_{0},\sigma_{0},p_{0})\in\mathbb{V}\mid
0\leq\varphi_{0}\leq 1,\,0\leq\sigma_{0}\leq\sigma_{\text{max}},\,0\leq
p_{0}\leq p_{\text{max}}\right\\},$ (2.19)

where $\sigma_{\text{max}},p_{\text{max}}\in L^{\infty}(\Omega)$ are given.
Note that $\mathcal{I}_{\text{ad}}$ is a closed and convex subset of
$\mathbb{V}\cap L^{\infty}(\Omega)^{3}$. Next, we want to study in detail the
operator $\mathcal{R}:\mathcal{I}_{\text{ad}}\to\mathbb{H}$ which associates
to any initial data $(\varphi_{0},\sigma_{0},p_{0})\in\mathcal{I}_{\text{ad}}$
the corresponding solution $(\varphi(T),\sigma(T),p(T))$ to (2.1)–(2.5),
evaluated at the final time, that is

$\mathcal{R}:\mathcal{I}_{\text{ad}}\to\mathbb{H},\quad\mathcal{R}((\varphi_{0},\sigma_{0},p_{0}))=(\varphi(T),\sigma(T),p(T)).$
(2.20)

Note that $\mathcal{R}$ is well-defined and Lipschitz-continuous by Theorem
2.1. We now want to prove that it is also continuously Fréchet-differentiable.
To do this, we introduce the linearised system:

$\displaystyle\partial_{t}Y-\lambda\Delta
Y+F^{\prime\prime}(\overline{\varphi})Y-m(\overline{\sigma})\mathbbm{h}^{\prime\prime}(\overline{\varphi})Y-m^{\prime}(\overline{\sigma})\mathbbm{h}^{\prime}(\overline{\varphi})Z=0\qquad$
in $Q_{T}$, (2.21) $\displaystyle\partial_{t}Z-\eta\Delta
Z+\gamma_{h}Z+\gamma_{ch}\overline{\sigma}Y+\gamma_{ch}\overline{\varphi}Z-S_{ch}Y=0$
in $Q_{T}$, (2.22) $\displaystyle\partial_{t}P-D\Delta
P+\gamma_{p}P=\alpha_{ch}Y$ in $Q_{T}$, (2.23) $\displaystyle
Y=0,\quad\partial_{\boldsymbol{n}}Z=\partial_{\boldsymbol{n}}P=0$ on
$\Sigma_{T}$, (2.24) $\displaystyle Y(0)=h,\quad Z(0)=k,\quad P(0)=w$ in
$\Omega$, (2.25)

where $h,k,w\in L^{2}(\Omega)$. This system is deduced by fixing an initial
datum $(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})$, as
well as its corresponding solution
$(\overline{\varphi},\overline{\sigma},\overline{p})$, by linearising near it,
i.e. by considering

$\varphi=\overline{\varphi}+Y,\,\sigma=\overline{\sigma}+Z,\,p=\overline{p}+P,\,\varphi_{0}=\overline{\varphi}_{0}+h,\,\sigma_{0}=\overline{\sigma}_{0}+k,\,p_{0}=\overline{p}_{0}+w,$

and by recovering the system solved by $(Y,Z,P)$, while approximating the non-
linearities at the first-order of their Taylor expansion. We have the
following result about the well-posedness of (2.21)–(2.25).

######  Proposition 2.4.

Assume hypotheses A1–A3 and let $(h,k,w)\in\mathbb{H}$. Let
$(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})\in\mathcal{I}_{\text{ad}}$.
Then, (2.21)–(2.25) admits a unique weak solution
$(Y,Z,P)\in\mathbb{Y}_{0}\times\mathbb{Y}\times\mathbb{Y}$, which solves the
system in variational formulation and satisfies the estimate:

$\lVert Y\rVert^{2}_{\mathbb{Y}_{0}}+\lVert Z\rVert^{2}_{\mathbb{Y}}+\lVert
P\rVert^{2}_{\mathbb{Y}}\leq C\left(\lVert h\rVert^{2}_{H}+\lVert
k\rVert^{2}_{H}+\lVert w\rVert^{2}_{H}\right).$

Moreover, if $(h,k,w)\in\mathbb{V}$, the solution enjoys higher regularity
$(Y,Z,P)\in\mathbb{X}_{0}\times\mathbb{X}\times\mathbb{X}$, solves the system
almost everywhere and satisfies the estimate:

$\lVert Y\rVert^{2}_{\mathbb{X}_{0}}+\lVert Z\rVert^{2}_{\mathbb{X}}+\lVert
P\rVert^{2}_{\mathbb{X}}\leq C\left(\lVert h\rVert^{2}_{V_{0}}+\lVert
k\rVert^{2}_{V}+\lVert w\rVert^{2}_{V}\right).$

###### Proof.

Note that system (2.21)–(2.25) is composed of three _linear_ parabolic
equations without sources and with coefficients in $L^{\infty}(Q_{T})$, thanks
to Theorem 2.1. Then, the statements of Proposition 2.4 follow by standard
results from the theory of parabolic equations. ∎

The linearised system (2.21)–(2.25) is an _ansatz_ for the explicit expression
of the Fréchet-derivative of the operator $\mathcal{R}$. Indeed, we are now in
position to show that

$\mathrm{D}\mathcal{R}(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})[(h,k,w)]=(Y(T),Z(T),P(T)).$
(2.26)

Let $\mathcal{I}_{R}$ be an open subset of $\mathbb{V}\cap
L^{\infty}(\Omega)^{3}$ such that
$\mathcal{I}_{\text{ad}}\subseteq\mathcal{I}_{R}$. Note that, in
$\mathcal{I}_{R}$, the continuous dependence estimate (2.12) holds with $C$
depending only on $R$ and the fixed data of the system. Our aim is to show
that $\mathcal{R}:\mathcal{I}_{R}\to\mathbb{H}$ is continuously Fréchet-
differentiable. Indeed, we can prove:

######  Theorem 2.5.

Assume hypothesis A1–A3. Then $\mathcal{R}:\mathcal{I}_{R}\to\mathbb{H}$ is
Fréchet-differentiable, i.e. for any
$(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})\in\mathcal{I}_{R}$
there exists a unique Fréchet-derivative
$\mathrm{D}\mathcal{R}((\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0}))\in\mathcal{L}(\mathbb{V}\cap
L^{\infty}(\Omega)^{3},\mathbb{H})$ such that, as
$\lVert(h,k,w)\rVert_{\mathbb{V}}\to 0$,

$\frac{\lVert\mathcal{R}((\overline{\varphi}_{0}+h,\overline{\sigma}_{0}+k,\overline{p}_{0}+w))-\mathcal{R}((\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0}))-\mathrm{D}\mathcal{R}((\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0}))[(h,k,w)]\rVert_{\mathbb{H}}}{\lVert(h,k,w)\rVert_{\mathbb{V}}}\to
0.$ (2.27)

Moreover, for any $(h,k,w)\in\mathbb{V}\cap L^{\infty}(\Omega)^{3}$, the
Fréchet-derivative at
$(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})$ in $(h,k,w)$
is defined as in (2.26), where $(Y(T),Z(T),P(T))$ is the solution to the
linearised system (2.21)–(2.25) with initial data $(h,k,w)$, evaluated at the
final time.

Additionally, the Fréchet-derivative $\mathrm{D}\mathcal{R}$ is Lipschitz-
continuous as a function from $\mathcal{I}_{R}$ to the space
$\mathcal{L}(\mathbb{V}\cap L^{\infty}(\Omega)^{3},\mathbb{H})$. More
precisely, the following estimate holds:

$\lVert\mathrm{D}\mathcal{R}((\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1}))-\mathrm{D}\mathcal{R}((\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2}))\rVert^{2}_{\mathcal{L}(\mathbb{V}\cap
L^{\infty}(\Omega)^{3},\mathbb{H})}\leq
C_{0}\lVert(\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1})-(\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2})\rVert^{2}_{\mathbb{H}},$
(2.28)

for some constant $C_{0}>0$ depending only on the parameters of the system.

###### Proof.

We observe that it is sufficient to prove the result for any small enough
perturbation $(h,k,w)$, i.e., we fix $\Lambda>0$ and consider only
perturbations such that

$\lVert(h,k,w)\rVert_{L^{\infty}(\Omega)^{3}}\leq\Lambda.$ (2.29)

Now, we fix $(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})$
and $(h,k,w)$ as above and consider

$\displaystyle(\varphi,\sigma,p)\text{ as the solution to
\eqref{eq:phi}--\eqref{ic} with data
}(\overline{\varphi}_{0}+h,\overline{\sigma}_{0}+k,\overline{p}_{0}+w),$
$\displaystyle(\overline{\varphi},\overline{\sigma},\overline{p})\text{ as the
solution to \eqref{eq:phi}--\eqref{ic} with data
}(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0}),$
$\displaystyle(Y,Z,P)\text{ as the solution to
\eqref{eq:philin}--\eqref{iclin} with data }(h,k,w).$

For the Fréchet-differentiability, it is enough to show that there exist a
constant $C>0$, depending only on the parameters of the system and possibly on
$\Lambda$, and an exponent $s>2$ such that

$\lVert(\varphi(T),\sigma(T),p(T))-(\overline{\varphi}(T),\overline{\sigma}(T),\overline{p}(T))-(Y(T),Z(T),P(T))\rVert^{2}_{\mathbb{H}}\leq
C\lVert(h,k,w)\rVert^{s}_{\mathbb{H}}.$

To do this, we introduce the additional variables

$\displaystyle\psi:=\varphi-\overline{\varphi}-Y\in\mathbb{X}_{0},$
$\displaystyle\theta:=\sigma-\overline{\sigma}-Z\in\mathbb{X},$
$\displaystyle\zeta:=p-\overline{p}-P\in\mathbb{X},$

which, by Theorem 2.1 and Proposition 2.4, enjoy the regularities shown above.
Then, this is equivalent to showing that

$\lVert(\psi(T),\theta(T),\zeta(T))\rVert_{\mathbb{H}}^{2}\leq
C\lVert(h,k,w)\rVert^{s}_{\mathbb{H}}.$ (2.30)

By inserting the equations solved by the variables in the definitions of
$\psi$, $\theta$ and $\zeta$ and exploiting the linearity of the involved
differential operators, we infer that $\psi$, $\theta$ and $\zeta$ satisfy:

$\displaystyle\partial_{t}\psi-\lambda\Delta\psi+F^{h}=M^{h}$
$\displaystyle\text{in }Q_{T},$ (2.31)
$\displaystyle\partial_{t}\theta-\eta\Delta\theta+\gamma_{h}\theta=Q^{h}+S_{ch}\psi\qquad$
$\displaystyle\text{in }Q_{T},$ (2.32)
$\displaystyle\partial_{t}\zeta-D\Delta\zeta+\gamma_{p}\zeta=\alpha_{ch}\psi$
$\displaystyle\text{in }Q_{T},$ (2.33)

together with boundary and initial conditions:

$\displaystyle\psi=0,\quad\partial_{\boldsymbol{n}}\theta=\partial_{\boldsymbol{n}}\zeta=0\qquad$
$\displaystyle\text{on }\Sigma_{T},$ (2.34)
$\displaystyle\psi(0)=0,\quad\theta(0)=0,\quad\zeta(0)=0\qquad$
$\displaystyle\text{in }\Omega,$ (2.35)

where

$\displaystyle F^{h}$
$\displaystyle=F^{\prime}(\varphi)-F^{\prime}(\overline{\varphi})-F^{\prime\prime}(\overline{\varphi})Y,$
$\displaystyle M^{h}$
$\displaystyle=m(\sigma)\mathbbm{h}^{\prime}(\varphi)-m(\overline{\sigma})\mathbbm{h}^{\prime}(\overline{\varphi})-m(\overline{\sigma})\mathbbm{h}^{\prime\prime}(\overline{\varphi})Y-m^{\prime}(\overline{\sigma})\mathbbm{h}^{\prime}(\overline{\varphi})Z,$
$\displaystyle Q^{h}$
$\displaystyle=\gamma_{ch}\sigma\varphi-\gamma_{ch}\overline{\sigma}\overline{\varphi}-\gamma_{ch}\overline{\sigma}Y-\gamma_{ch}\overline{\varphi}Z.$

Before going on, we can rewrite in a better way the terms $F^{h}$, $M^{h}$ and
$Q^{h}$, by using the following version of Taylor’s theorem with integral
remainder for a real function $f\in\mathcal{C}^{2}$ at a point
$x_{0}\in\mathbb{R}$:

$f(x)=f(x_{0})+f^{\prime}(x_{0})(x-x_{0})+\left(\int_{0}^{1}(1-z)f^{\prime\prime}(x_{0}+z(x-x_{0}))\,\mathrm{d}z\right)(x-x_{0})^{2}.$

Indeed, with straightforward calculations, one can see that

$\displaystyle
F^{h}=F^{\prime\prime}(\overline{\varphi})\psi+R_{1}^{h}(\varphi-\overline{\varphi})^{2},$

and also, up to adding and subtracting some additional terms, that

$\displaystyle M^{h}$
$\displaystyle=(m(\sigma)-m(\overline{\sigma}))(\mathbbm{h}^{\prime}(\varphi)-\mathbbm{h}^{\prime}(\overline{\varphi}))+m(\overline{\sigma})[\mathbbm{h}^{\prime\prime}(\overline{\varphi})\psi+R_{2}^{h}(\varphi-\overline{\varphi})^{2}]$
$\displaystyle\quad+\mathbbm{h}^{\prime}(\varphi)[m^{\prime}(\overline{\sigma})\theta+R_{3}^{h}(\sigma-\overline{\sigma})^{2}],$
$\displaystyle Q^{h}$
$\displaystyle=\gamma_{ch}(\sigma-\overline{\sigma})(\varphi-\overline{\varphi})+\gamma_{ch}\overline{\sigma}\psi+\gamma_{ch}\overline{\varphi}\theta,$

where

$\displaystyle
R_{1}^{h}=\int_{0}^{1}(1-z)F^{\prime\prime\prime}(\overline{\varphi}+z(\varphi-\overline{\varphi}))\,\mathrm{d}z,\quad
R_{2}^{h}=\int_{0}^{1}(1-z)\mathbbm{h}^{\prime\prime\prime}(\overline{\varphi}+z(\varphi-\overline{\varphi}))\,\mathrm{d}z,$
$\displaystyle
R_{3}^{h}=\int_{0}^{1}(1-z)m^{\prime\prime}(\overline{\sigma}+z(\sigma-\overline{\sigma}))\,\mathrm{d}z.$

By using global boundedness of the variables, given by Theorem 2.1, one can
easily see that there exists a constant $C>0$, depending only on the
parameters of the system and on $\Lambda$, such that

$\lVert R_{1}^{h}\rVert_{L^{\infty}(Q_{T})},\lVert
R_{2}^{h}\rVert_{L^{\infty}(Q_{T})},\lVert
R_{3}^{h}\rVert_{L^{\infty}(Q_{T})}\leq C.$ (2.36)

In order to show (2.30), we now proceed with a priori estimates on the system
(2.31)–(2.35). Indeed, we test (2.31) by $\psi$, (2.32) by $\theta$, (2.33) by
$\zeta$ and sum them up to obtain:

$\begin{split}&\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\lVert\varphi\rVert^{2}_{H}+\lVert\sigma\rVert^{2}_{H}+\lVert
p\rVert^{2}_{H}\right)+\lambda\lVert\nabla\varphi\rVert^{2}_{H}+\eta\lVert\nabla\sigma\rVert^{2}_{H}+D\lVert\nabla
p\rVert^{2}_{H}+\gamma_{h}\lVert\theta\rVert^{2}_{H}+\gamma_{p}\lVert\zeta\rVert^{2}_{H}\\\
&\quad=-(F^{h},\psi)_{H}+(M^{h},\psi)_{H}+(Q^{h},\theta)_{H}+S_{ch}(\psi,\theta)_{H}+\alpha_{ch}(\psi,\zeta)_{H}.\end{split}$
(2.37)

The last two terms on the right-hand side of (2.37) can easily be treated by
means of Cauchy-Schwarz and Young’s inequalities, therefore we focus on
estimating the first three. Indeed, by using again Hölder and Young’s
inequalities, together with (2.36), the global boundedness of
$(\overline{\varphi},\overline{\sigma},\overline{p})$ given by Theorem 2.1,
the regularity and local Lipschitz continuity of $F$, $m$ and $\mathbbm{h}$
and the Sobolev embedding $V\hookrightarrow L^{4}(\Omega)$, we infer that

$\displaystyle(F^{h},\psi)_{H}$
$\displaystyle=\int_{\Omega}F^{\prime\prime}(\overline{\varphi})\psi^{2}+R_{1}^{h}(\varphi-\overline{\varphi})^{2}\psi\,\mathrm{d}x\leq
C\lVert\psi\rVert^{2}_{H}+C\lVert\varphi-\overline{\varphi}\rVert^{4}_{L^{4}(\Omega)}$
$\displaystyle\leq
C\lVert\psi\rVert^{2}_{H}+C\lVert\varphi-\overline{\varphi}\rVert^{4}_{V_{0}},$
$\displaystyle(M^{h},\psi)_{H}$
$\displaystyle=\int_{\Omega}(m(\sigma)-m(\overline{\sigma}))(\mathbbm{h}^{\prime}(\varphi)-\mathbbm{h}^{\prime}(\overline{\varphi}))\psi\,\mathrm{d}x+\int_{\Omega}m(\overline{\sigma})[\mathbbm{h}^{\prime\prime}(\overline{\varphi})\psi+R_{2}^{h}(\varphi-\overline{\varphi})^{2}]\psi\,\mathrm{d}x$
$\displaystyle\quad+\int_{\Omega}\mathbbm{h}^{\prime}(\varphi)[m^{\prime}(\overline{\sigma})\theta+R_{3}^{h}(\sigma-\overline{\sigma})^{2}]\psi\,\mathrm{d}x$
$\displaystyle\leq
C\lVert\psi\rVert^{2}_{H}+C\lVert\theta\rVert^{2}_{H}+C\lVert\varphi-\overline{\varphi}\rVert^{4}_{L^{4}(\Omega)}+C\lVert\sigma-\overline{\sigma}\rVert^{4}_{L^{4}(\Omega)}$
$\displaystyle\leq
C\lVert\psi\rVert^{2}_{H}+C\lVert\theta\rVert^{2}_{H}+C\lVert\varphi-\overline{\varphi}\rVert^{4}_{V_{0}}+C\lVert\sigma-\overline{\sigma}\rVert^{4}_{V},$
$\displaystyle(Q^{h},\theta)_{H}$
$\displaystyle=\gamma_{ch}\int_{\Omega}(\sigma-\overline{\sigma})(\varphi-\overline{\varphi})\theta+\overline{\sigma}\psi\theta+\overline{\varphi}\theta^{2}\,\mathrm{d}x$
$\displaystyle\leq
C\lVert\theta\rVert^{2}_{H}+C\lVert\psi\rVert^{2}_{H}+C\lVert\sigma-\overline{\sigma}\rVert^{2}_{L^{4}(\Omega)}\lVert\varphi-\overline{\varphi}\rVert^{2}_{L^{4}(\Omega)}$
$\displaystyle\leq
C\lVert\psi\rVert^{2}_{H}+C\lVert\theta\rVert^{2}_{H}+C\lVert\varphi-\overline{\varphi}\rVert^{4}_{V_{0}}+C\lVert\sigma-\overline{\sigma}\rVert^{4}_{V}.$

Then, by integrating on $(0,t)$ for any $t\in(0,T)$, we get that

$\displaystyle\frac{1}{2}\lVert\psi(t)\rVert^{2}_{H}+\frac{1}{2}\lVert\theta(t)\rVert^{2}_{H}+\frac{1}{2}\lVert\zeta(t)\rVert^{2}_{H}+\lambda\int_{0}^{t}\lVert\nabla\psi\rVert^{2}_{H}\,\mathrm{d}s+\eta\int_{0}^{t}\lVert\nabla\theta\rVert^{2}_{H}\,\mathrm{d}s+D\int_{0}^{t}\lVert\nabla\zeta\rVert^{2}_{H}\,\mathrm{d}s$
$\displaystyle\quad\leq
C\int_{0}^{T}\lVert\psi\rVert^{2}_{H}\,\mathrm{d}s+C\int_{0}^{T}\lVert\theta\rVert^{2}_{H}\,\mathrm{d}s+C\int_{0}^{T}\lVert\zeta\rVert^{2}_{H}\,\mathrm{d}s+C\int_{0}^{T}\lVert\varphi-\overline{\varphi}\rVert^{4}_{V_{0}}\,\mathrm{d}s+C\int_{0}^{T}\lVert\sigma-\overline{\sigma}\rVert^{4}_{V}\,\mathrm{d}s.$

Therefore, by applying Gronwall’s lemma and the continuous dependence estimate
found in Proposition 2.3, we conclude that

$\begin{split}&\lVert\psi\rVert^{2}_{L^{\infty}(0,T;H)\cap
L^{2}(0,T;V_{0})}+\lVert\theta\rVert^{2}_{L^{\infty}(0,T;H)\cap
L^{2}(0,T;V)}\\\ &\quad+\lVert\zeta\rVert^{2}_{L^{\infty}(0,T;H)\cap
L^{2}(0,T;V)}\leq C\left(\lVert h\rVert^{4}_{V_{0}}+\lVert
k\rVert^{4}_{V}+\lVert w\rVert^{4}_{V}\right).\end{split}$ (2.38)

Moreover, by comparison in (2.31), (2.32) and (2.33), starting from (2.38), we
also deduce that

$\lVert\psi\rVert^{2}_{H^{1}(0,T;V_{0}^{*})}+\lVert\theta\rVert^{2}_{H^{1}(0,T;V^{*})}+\lVert\zeta\rVert^{2}_{H^{1}(0,T;V^{*})}\leq
C\left(\lVert h\rVert^{4}_{V_{0}}+\lVert k\rVert^{4}_{V}+\lVert
w\rVert^{4}_{V}\right).$ (2.39)

Therefore, due to the standard embedding $H^{1}(0,T;\mathbb{V}^{*})\cap
L^{2}(0,T;\mathbb{V})\hookrightarrow\mathcal{C}^{0}([0,T];\mathbb{H})$,
estimates (2.38) and (2.39) imply that

$\lVert(\psi(T),\theta(T),\zeta(T))\rVert^{2}_{\mathbb{H}}\leq
C\lVert(h,k,w)\rVert^{4}_{\mathbb{V}},$

which is exactly (2.30) with $s=4>2$. Then, the operator $\mathcal{R}$ is
Fréchet-differentiable and its derivative is fully characterised.

Next, we want to prove that the Fréchet derivative is Lipschitz continuous. To
do this, we show that, given two sets of initial data
$(\overline{\varphi}_{0}^{i},\overline{\sigma}_{0}^{i},\overline{p}_{0}^{i})\in\mathcal{I}_{\text{ad}}$,
$i=1,2$, with corresponding solutions
$(\overline{\varphi}_{i},\overline{\sigma}_{i},\overline{p}_{i})$, $i=1,2$,
then for any $(h,k,w)\in\mathbb{V}\cap L^{\infty}(\Omega)^{3}$ it holds that

$\begin{split}&\lVert\mathrm{D}\mathcal{R}((\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1}))[(h,k,w)]-\mathrm{D}\mathcal{R}((\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2}))[(h,k,w)]\rVert^{2}_{\mathbb{H}}\\\
&\quad=\lVert(Y_{1}(T),Z_{1}(T),P_{1}(T))-(Y_{2}(T),Z_{2}(T),P_{2}(T))\rVert^{2}_{\mathbb{H}}\\\
&\leq
C_{0}\lVert(\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1})-(\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2})\rVert^{2}_{\mathbb{H}}.\end{split}$
(2.40)

To prove (2.40), we consider the system solved by the differences
$Y=Y_{1}-Y_{2}$, $Z=Z_{1}-Z_{2}$ and $P=P_{1}-P_{2}$ of the corresponding
solutions of the linearised system (2.21)–(2.25), which, up to adding and
subtracting some terms, has the following form:

$\displaystyle\partial_{t}Y-\lambda\Delta
Y+F^{\prime\prime}(\overline{\varphi}_{1})Y+(F^{\prime\prime}(\overline{\varphi}_{1})-F^{\prime\prime}(\overline{\varphi}_{2}))Y_{2}$
$\displaystyle\quad=m(\overline{\sigma}_{1})\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y+(m(\overline{\sigma}_{1})-m(\overline{\sigma}_{2}))\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y_{2}$
$\displaystyle\qquad+m(\overline{\sigma}_{2})(\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{2}))Y_{2}+m^{\prime}(\overline{\sigma}_{1})\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z$
$\displaystyle\qquad+(m^{\prime}(\overline{\sigma}_{1})-m^{\prime}(\overline{\sigma}_{2}))\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z_{2}+m^{\prime}(\overline{\sigma}_{1})(\mathbbm{h}^{\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime}(\overline{\varphi}_{2}))Z_{2}\qquad$
in $Q_{T}$, (2.41) $\displaystyle\partial_{t}Z-\eta\Delta
Z+\gamma_{h}Z+\gamma_{ch}\overline{\sigma}_{1}Y+\gamma_{ch}(\overline{\sigma}_{1}-\overline{\sigma}_{2})Y_{2}$
$\displaystyle\qquad+\gamma_{ch}\overline{\varphi}_{1}Z+\gamma_{ch}(\overline{\varphi}_{1}-\overline{\varphi}_{2})Z_{2}=S_{ch}Y$
in $Q_{T}$, (2.42) $\displaystyle\partial_{t}P-D\Delta
P+\gamma_{p}P=\alpha_{ch}Y$ in $Q_{T}$, (2.43) $\displaystyle
Y=0,\quad\partial_{\boldsymbol{n}}Z=\partial_{\boldsymbol{n}}P=0$ on
$\Sigma_{T}$, (2.44) $\displaystyle Y(0)=0,\quad Z(0)=0,\quad P(0)=0$ in
$\Omega$. (2.45)

Then, we test (2.41) by $Y$, (2.42) by $Z$, (2.43) by $P$ and sum them up to
obtain:

$\begin{split}&\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\lVert
Y\rVert^{2}_{H}+\lVert Z\rVert^{2}_{H}+\lVert P\rVert^{2}_{H}\right)\\\
&\qquad+\lambda\lVert\nabla Y\rVert^{2}_{H}+\eta\lVert\nabla
Z\rVert^{2}_{H}+D\lVert\nabla P\rVert^{2}_{H}+\gamma_{h}\lVert
Z\rVert^{2}_{H}+\gamma_{p}\lVert P\rVert^{2}_{H}\\\
&\quad=(F^{\prime\prime}(\overline{\varphi}_{1})Y,Y)_{H}+((F^{\prime\prime}(\overline{\varphi}_{1})-F^{\prime\prime}(\overline{\varphi}_{2}))Y_{2},Y)_{H}+(m(\overline{\sigma}_{1})\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y,Y)_{H}\\\
&\qquad+((m(\overline{\sigma}_{1})-m(\overline{\sigma}_{2}))\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y_{2},Y)_{H}+(m(\overline{\sigma}_{2})(\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{2}))Y_{2},Y)_{H}\\\
&\qquad+(m^{\prime}(\overline{\sigma}_{1})\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z,Y)+((m^{\prime}(\overline{\sigma}_{1})-m^{\prime}(\overline{\sigma}_{2}))\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z_{2},Y)_{H}\\\
&\qquad+(m^{\prime}(\overline{\sigma}_{1})(\mathbbm{h}^{\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime}(\overline{\varphi}_{2}))Z_{2},Y)_{H}-\gamma_{ch}(\overline{\sigma}_{1}Y,Z)_{H}-\gamma_{ch}((\overline{\sigma}_{1}-\overline{\sigma}_{2})Y_{2},Z)_{H}\\\
&\qquad-\gamma_{ch}(\overline{\varphi}_{1}Z,Z)_{H}-\gamma_{ch}((\overline{\varphi}_{1}-\overline{\varphi}_{2})Z_{2},Z)_{H}+(S_{ch}Y,Z)_{H}+\alpha_{ch}(Y,P)_{H}.\end{split}$
(2.46)

Now, by using Hölder and Young’s inequalities, together with the global
boundedness of
$(\overline{\varphi}_{i},\overline{\sigma}_{i},\overline{p}_{i})$ given by
Theorem 2.1, the regularity and local Lipschitz continuity of $F$, $m$ and
$\mathbbm{h}$ and the Sobolev embedding $V\hookrightarrow L^{4}(\Omega)$, we
can easily estimate all the terms on the right-hand side of (2.46) in the
following way:

$\displaystyle(F^{\prime\prime}(\overline{\varphi}_{1})Y,Y)_{H}+(m(\overline{\sigma}_{1})\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y,Y)_{H}+(m^{\prime}(\overline{\sigma}_{1})\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z,Y)$
$\displaystyle\qquad-\gamma_{ch}(\overline{\sigma}_{1}Y,Z)_{H}-\gamma_{ch}(\overline{\varphi}_{1}Z,Z)_{H}+(S_{ch}Y,Z)_{H}+\alpha_{ch}(Y,P)_{H}$
$\displaystyle\quad\leq C\lVert Y\rVert^{2}_{H}+C\lVert
Z\rVert^{2}_{H}+C\lVert P\rVert^{2}_{H},$
$\displaystyle((F^{\prime\prime}(\overline{\varphi}_{1})-F^{\prime\prime}(\overline{\varphi}_{2}))Y_{2},Y)_{H}+((m(\overline{\sigma}_{1})-m(\overline{\sigma}_{2}))\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})Y_{2},Y)_{H}$
$\displaystyle\qquad+(m(\overline{\sigma}_{2})(\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime\prime}(\overline{\varphi}_{2}))Y_{2},Y)_{H}+((m^{\prime}(\overline{\sigma}_{1})-m^{\prime}(\overline{\sigma}_{2}))\mathbbm{h}^{\prime}(\overline{\varphi}_{1})Z_{2},Y)_{H}$
$\displaystyle\qquad+(m^{\prime}(\overline{\sigma}_{1})(\mathbbm{h}^{\prime}(\overline{\varphi}_{1})-\mathbbm{h}^{\prime}(\overline{\varphi}_{2}))Z_{2},Y)_{H}-\gamma_{ch}((\overline{\sigma}_{1}-\overline{\sigma}_{2})Y_{2},Z)_{H}-\gamma_{ch}((\overline{\varphi}_{1}-\overline{\varphi}_{2})Z_{2},Z)_{H}$
$\displaystyle\quad\leq
C\lVert\overline{\varphi}_{1}-\overline{\varphi}_{2}\rVert_{L^{4}(\Omega)}\lVert
Y_{2}\rVert_{H}\lVert
Y\rVert_{L^{4}(\Omega)}+C\lVert\overline{\sigma}_{1}-\overline{\sigma}_{2}\rVert_{L^{4}(\Omega)}\lVert
Y_{2}\rVert_{H}\lVert Y\rVert_{L^{4}(\Omega)}$
$\displaystyle\qquad+C\lVert\overline{\sigma}_{1}-\overline{\sigma}_{2}\rVert_{L^{4}(\Omega)}\lVert
Z_{2}\rVert_{H}\lVert
Y\rVert_{L^{4}(\Omega)}+C\lVert\overline{\varphi}_{1}-\overline{\varphi}_{2}\rVert_{L^{4}(\Omega)}\lVert
Z_{2}\rVert_{H}\lVert Y\rVert_{L^{4}(\Omega)}$
$\displaystyle\qquad+C\lVert\overline{\sigma}_{1}-\overline{\sigma}_{2}\rVert_{L^{4}(\Omega)}\lVert
Y_{2}\rVert_{H}\lVert
Z\rVert_{L^{4}(\Omega)}+C\lVert\overline{\varphi}_{1}-\overline{\varphi}_{2}\rVert_{L^{4}(\Omega)}\lVert
Z_{2}\rVert_{H}\lVert Z\rVert_{L^{4}(\Omega)}$
$\displaystyle\quad\leq\frac{\lambda}{2}\lVert
Y\rVert^{2}_{V_{0}}+\frac{\eta}{2}\lVert
Z\rVert^{2}_{V}+C\underbrace{\left(\lVert Y_{2}\rVert^{2}_{H}+\lVert
Z_{2}\rVert^{2}_{H}\right)}_{\in\,L^{\infty}(0,T)}\lVert\overline{\varphi}_{1}-\overline{\varphi}_{2}\rVert^{2}_{V_{0}}+C\underbrace{\left(\lVert
Y_{2}\rVert^{2}_{H}+\lVert
Z_{2}\rVert^{2}_{H}\right)}_{\in\,L^{\infty}(0,T)}\lVert\overline{\sigma}_{1}-\overline{\sigma}_{2}\rVert^{2}_{V},$

where $\lVert Y_{2}\rVert^{2}_{H}$ and $\lVert Z_{2}\rVert^{2}_{H}$ are
uniformly bounded in $L^{\infty}(0,T)$ by Proposition 2.4. Therefore, by also
integrating on $(0,t)$, for any $t\in(0,T)$, starting from (2.46), we deduce
that

$\displaystyle\frac{1}{2}\left(\lVert Y(t)\rVert^{2}_{H}+\lVert
Z(t)\rVert^{2}_{H}+\lVert
P(t)\rVert^{2}_{H}\right)+\frac{\lambda}{2}\int_{0}^{t}\lVert\nabla
Y\rVert^{2}_{H}\,\mathrm{d}s+\frac{\eta}{2}\int_{0}^{t}\lVert\nabla
Z\rVert^{2}_{H}\,\mathrm{d}s+\frac{D}{2}\int_{0}^{t}\lVert\nabla
P\rVert^{2}_{H}\,\mathrm{d}s$ $\displaystyle\quad\leq C\int_{0}^{T}\lVert
Y\rVert^{2}_{H}+\lVert Z\rVert^{2}_{H}+\lVert
P\rVert^{2}_{H}\,\mathrm{d}s+C\int_{0}^{T}\lVert\overline{\varphi}_{1}-\overline{\varphi}_{2}\rVert^{2}_{V_{0}}+\lVert\overline{\sigma}_{1}-\overline{\sigma}_{2}\rVert^{2}_{V}\,\mathrm{d}s.$

Hence, by using Gronwall’s lemma and the continuous dependence estimate
(2.10), we infer that

$\begin{split}&\lVert Y\rVert^{2}_{L^{\infty}(0,T;H)\cap
L^{2}(0,T;V_{0})}+\lVert Z\rVert^{2}_{L^{\infty}(0,T;H)\cap L^{2}(0,T;V)}\\\
&\quad+\lVert P\rVert^{2}_{L^{\infty}(0,T;H)\cap L^{2}(0,T;V)}\leq
C\lVert(\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1})-(\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2})\rVert^{2}_{\mathbb{H}}.\end{split}$
(2.47)

Moreover, by comparison in (2.41), (2.42) and (2.43), starting from (2.47), we
also deduce that

$\lVert Y\rVert^{2}_{H^{1}(0,T;V_{0}^{*})}+\lVert
Z\rVert^{2}_{H^{1}(0,T;V^{*})}+\lVert P\rVert^{2}_{H^{1}(0,T;V^{*})}\leq
C\lVert(\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1})-(\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2})\rVert^{2}_{\mathbb{H}}.$
(2.48)

Therefore, due to the standard embedding $H^{1}(0,T;\mathbb{V}^{*})\cap
L^{2}(0,T;\mathbb{V})\hookrightarrow\mathcal{C}([0,T];\mathbb{H})$, estimates
(2.47) and (2.48) allow us to say that

$\lVert(Y(T),Z(T),P(T))\rVert^{2}_{\mathbb{H}}\leq
C\lVert(\overline{\varphi}_{0}^{1},\overline{\sigma}_{0}^{1},\overline{p}_{0}^{1})-(\overline{\varphi}_{0}^{2},\overline{\sigma}_{0}^{2},\overline{p}_{0}^{2})\rVert^{2}_{\mathbb{H}},$

which immediately implies (2.40). Note that all this procedure is independent
of the particular choice of $(h,k,w)\in\mathbb{V}\cap L^{\infty}(\Omega)^{3}$,
thus by taking the supremum in $(h,k,w)$ one recovers exactly (2.28). This
concludes the proof of Theorem 2.5. ∎

######  Remark 2.6.

Note that, by Proposition 2.4, the solution to the linearised system
(2.21)–(2.25) is well-defined for any $(h,k,w)\in\mathbb{V}$, without
additionally asking that $(h,k,w)\in\mathbb{V}\cap L^{\infty}(\Omega)^{3}$.
This means that, for any fixed
$(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})\in\mathcal{I}_{\text{ad}}$,
the Fréchet-derivative
$\mathrm{D}\mathcal{R}(\overline{\varphi}_{0},\overline{\sigma}_{0},\overline{p}_{0})$
can be extended to a continuous linear operator on the whole $\mathbb{V}$.

##  3 Analysis of the inverse problem

We devote this section to the analysis of the inverse problem of
reconstructing the initial data, given a measurement at the final time. The
first step would be to prove a backward uniqueness result for the system
(2.1)–(2.5), namely to show that if $(\varphi_{1},\sigma_{1},p_{1})$ and
$(\varphi_{1},\sigma_{1},p_{1})$ are two solutions such that
$(\varphi_{1},\sigma_{1},p_{1})(T)=(\varphi_{2},\sigma_{2},p_{2})(T)$, then
$(\varphi_{1},\sigma_{1},p_{1})(t)=(\varphi_{2},\sigma_{2},p_{2})(t)$ for any
$t\in[0,T]$. In particular, this would imply that the map $\mathcal{R}$ is
injective i.e. uniqueness for the inverse problem in
$\mathcal{I}_{\text{ad}}$. However, following a logarithmic convexity
approach, we prove a conditional stability result which consequently implies
uniqueness as a byproduct. Namely, assuming the initial data to belong to
$\mathcal{I}_{\text{ad}}$, we prove an Hölder stability estimate for
$\lVert(\varphi_{1}(t),\sigma_{1}(t),p_{1}(t))-(\varphi_{2}(t),\sigma_{2}(t),p_{2}(t)\rVert_{\mathbb{H}}$
for any $t>0$. To establish this, we exploit the global bounds on the
solutions of (2.1)–(2.5), obtained in Theorem 2.1. In particular, in the
sequel we will denote by $M>0$ the minimal constant such that

$\lVert(\varphi,\sigma,p)\rVert_{\mathcal{C}^{0}([0,T];\mathbb{H})}\leq M,$
(3.1)

uniformly for $(\varphi_{0},\sigma_{0},p_{0})\in\mathcal{I}_{\text{ad}}$. To
prove our conditional stability result, we follow [33, Theorem 3.1.3], which
is actually a simpler version of the results contained in [1]. Indeed, we can
prove:

######  Proposition 3.1.

Assume hypotheses A1–A3. Let $(\varphi_{1},\sigma_{1},p_{1})$ and
$(\varphi_{2},\sigma_{2},p_{2})$ be two solutions of (2.1)–(2.5) corresponding
to two triples of initial data
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})\in\mathcal{I}_{\text{ad}}$ for
$i=1,2$.

Then, there exists a constant $C_{1}>0$, depending only on the parameters of
the system, such that for any $t\in[0,T]$ the following estimate holds:

$\begin{split}&\lVert(\varphi_{1}(t),\sigma_{1}(t),p_{1}(t))-(\varphi_{2}(t),\sigma_{2}(t),p_{2}(t))\rVert_{\mathbb{H}}\\\
&\quad\leq
C_{1}M^{1-\lambda(t)}\lVert(\varphi_{1}(T),\sigma_{1}(T),p_{1}(T))-(\varphi_{2}(T),\sigma_{2}(T),p_{2}(T))\rVert_{\mathbb{H}}^{\lambda(t)}.\end{split}$
(3.2)

where $\lambda=\lambda(t)\geq 0$ can be either

$\lambda=\lambda_{1}(t)=\frac{1-e^{-\gamma t}}{1-e^{-\gamma
T}}\quad\text{or}\quad\lambda=\lambda_{2}(t)=\frac{e^{\gamma t}-1}{e^{\gamma
T}-1},$ (3.3)

for any $t\in[0,T]$, with $\gamma>0$ depending only on the fixed parameters of
the system.

###### Proof.

We call $\varphi=\varphi_{1}-\varphi_{2}$, $\sigma=\sigma_{1}-\sigma_{2}$,
$p=p_{1}-p_{2}$ and observe that they satisfy the system:

$\displaystyle\partial_{t}\varphi-\lambda\Delta\varphi=f_{\varphi}$ in
$Q_{T}$, (3.4) $\displaystyle\partial_{t}\sigma-\eta\Delta\sigma=f_{\sigma}$
in $Q_{T}$, (3.5) $\displaystyle\partial_{t}p-D\Delta p=f_{p}$ in $Q_{T}$,
(3.6)
$\displaystyle\varphi=0,\quad\partial_{\boldsymbol{n}}\sigma=\partial_{\boldsymbol{n}}p=0$
in $\Sigma_{T}$, (3.7)
$\displaystyle\varphi(0)=\varphi_{0}^{1}-\varphi_{0}^{2},\quad\sigma(0)=\sigma_{0}^{1}-\sigma_{0}^{2},\quad
p(0)=p_{0}^{1}-p_{0}^{2}\qquad$ in $\Omega$, (3.8)

where

$\displaystyle
f_{\varphi}:=-(F^{\prime}(\varphi_{1})-F^{\prime}(\varphi_{2}))+m(\sigma_{1})(\mathbbm{h}^{\prime}(\varphi_{1})-\mathbbm{h}^{\prime}(\varphi_{2}))+(m(\sigma_{1})-m(\sigma_{2}))\mathbbm{h}^{\prime}(\varphi_{2}),$
$\displaystyle
f_{\sigma}:=-\gamma_{h}\sigma-\gamma_{ch}\sigma_{1}(\varphi_{1}-\varphi_{2})-\gamma_{ch}(\sigma_{1}-\sigma_{2})\varphi_{2}+S_{ch}\varphi,$
$\displaystyle f_{p}:=-\gamma_{p}p+\alpha_{ch}\varphi.$

We start by estimating the terms on the right-hand side. Indeed, by using the
fact that $\varphi_{i},\sigma_{i}$, $i=1,2$, are uniformly bounded in
$L^{\infty}(Q_{T})$ and $F,\mathbbm{h}\in\mathcal{C}^{2}(\mathbb{R})$,
$m\in\mathcal{C}^{1}(\mathbb{R})$ and thus locally Lipschitz, we can infer
that

$\displaystyle\lVert f_{\varphi}\rVert_{H}$ $\displaystyle\leq\lVert
F^{\prime}(\varphi_{1})-F^{\prime}(\varphi_{2})\rVert_{H}+\lVert
m(\sigma_{1})\rVert_{\infty}\lVert\mathbbm{h}^{\prime}(\varphi_{1})-\mathbbm{h}^{\prime}(\varphi_{2})\rVert_{H}+\lVert\mathbbm{h}^{\prime}(\varphi_{2})\rVert_{\infty}\lVert
m(\sigma_{1})-m(\sigma_{2})\rVert_{H}$ $\displaystyle\leq
C\left(\lVert\varphi\rVert_{H}+\lVert\sigma\rVert_{H}\right).$

Similarly, we can also deduce that

$\displaystyle\lVert f_{\sigma}\rVert_{H}$
$\displaystyle\leq\gamma_{h}\lVert\sigma\rVert_{H}+\gamma_{ch}\lVert\sigma_{1}\rVert_{\infty}\lVert\varphi_{1}-\varphi_{2}\rVert_{H}+\gamma_{ch}\lVert\varphi_{2}\rVert_{\infty}\lVert\sigma_{1}-\sigma_{2}\rVert_{H}+S_{ch}\lVert\varphi\rVert_{H}$
$\displaystyle\leq
C\left(\lVert\sigma\rVert_{H}+\lVert\varphi\rVert_{H}\right),$

and that

$\lVert f_{p}\rVert_{H}\leq C\left(\lVert
p\rVert_{H}+\lVert\varphi\rVert_{H}\right).$

Therefore, in the end, we showed that

$\lVert f_{\varphi}\rVert_{H}+\lVert f_{\sigma}\rVert_{H}+\lVert
f_{p}\rVert_{H}\leq
C\left(\lVert\varphi\rVert_{H}+\lVert\sigma\rVert_{H}+\lVert
p\rVert_{H}\right),\quad\text{for a.e. }t\in[0,T].$ (3.9)

To simplify notations, we now introduce a vector-valued formulation by
defining:

$\displaystyle\boldsymbol{\uppsi}:=(\varphi,\sigma,p),\quad\nabla\boldsymbol{\uppsi}:=(\nabla\varphi,\nabla\sigma,\nabla
p),$
$\displaystyle\Delta\boldsymbol{\uppsi}:=(\Delta\varphi,\Delta\sigma,\Delta
p),\quad
f_{\boldsymbol{\uppsi}}:=(f_{\varphi},f_{\sigma},f_{p}),\quad\boldsymbol{D}:=\operatorname{diag}(\lambda,\eta,D).$

Then, the previous estimate (3.9) now becomes

$\lVert f_{\boldsymbol{\uppsi}}\rVert_{\mathbb{H}}\leq
C\lVert\boldsymbol{\uppsi}\rVert_{\mathbb{H}}.$ (3.10)

Thus, we have shown that the system solved by $(\varphi,\sigma,p)$ has the
form

$\partial_{t}\boldsymbol{\uppsi}-\boldsymbol{D}\Delta\boldsymbol{\uppsi}=f_{\boldsymbol{\uppsi}},\hbox{
with $\lVert f_{\boldsymbol{\uppsi}}\rVert_{\mathbb{H}}\leq
C\lVert\boldsymbol{\uppsi}\rVert_{\mathbb{H}}$},$

where $C$ depends only on the parameters of the system and on
$\lVert\varphi_{0}^{i}\rVert_{\infty},\lVert\sigma_{0}^{i}\rVert_{\infty}$,
$i=1,2$. Next, recall that
$\boldsymbol{\uppsi}\in\mathbb{X}_{0}\times\mathbb{X}\times\mathbb{X}$,
therefore, by respectively testing (3.4) by $\varphi$ and $-\Delta\varphi$,
(3.5) by $\sigma$ and $-\Delta\sigma$, (3.6) by $p$ and $-\Delta p$ and
summing them up, we obtain the identities:

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\lVert\boldsymbol{\uppsi}\rVert^{2}_{\mathbb{H}}+\lVert\boldsymbol{D}\nabla\boldsymbol{\uppsi}\rVert^{2}_{\mathbb{H}}=(f_{\boldsymbol{\uppsi}},\boldsymbol{\uppsi})_{\mathbb{H}},$
(3.11)
$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\lVert\nabla\boldsymbol{\uppsi}\rVert^{2}_{\mathbb{H}}+\lVert\boldsymbol{D}\Delta\boldsymbol{\uppsi}\rVert^{2}_{\mathbb{H}}=(f_{\boldsymbol{\uppsi}},-\Delta\boldsymbol{\uppsi})_{\mathbb{H}}.$
(3.12)

Now, we are exactly in the situation to apply [33, Theorem 3.1.3], with
$\alpha=C$, to obtain the desired estimate. In particular, observe that, by
going through the proof of the cited Theorem and especially the one of [33,
Lemma 3.1.4], the argument can be repeated even in our regularity setting with
$\boldsymbol{\uppsi}\in\mathbb{X}_{0}\times\mathbb{X}\times\mathbb{X}$.
Indeed, the crucial point is the validity of the identities (3.11) and (3.12).
For the sake of completeness, we recall here below the fundamental steps of
the proof.

Let

$q:=\lVert\boldsymbol{\uppsi}\rVert^{2}_{\mathbb{H}},\quad
g:=\frac{2(f_{\boldsymbol{\uppsi}},\boldsymbol{\uppsi})_{\mathbb{H}}}{q},\quad
l:=\log q-\int_{0}^{t}g\,\mathrm{d}s.$

Then, the first step is to show, through explicit calculations, that there
exists a constant $c\geq 0$ such that $l$ satisfies the differential
inequality

$\partial^{2}_{t}l+c\lvert\partial_{t}l\rvert+c\geq 0\quad\text{in $(0,T)$}.$
(3.13)

This can be done exactly as in [33, Lemma 3.1.4], by relying on the identities
(3.11) and (3.12). Note that these calculations are actually easier in our
case, since our operator $A=-\boldsymbol{D}\Delta$ has no skew-symmetric part.
Then, exactly as in [33, Lemma 3.1.5], by exploiting maximum, minimum and
comparison principles, we can see that $l(t)\leq L(t)$ for any $t\in(0,T)$,
where $L$ is the solution to the following boundary value problem:

$\begin{cases}\partial_{t}^{2}L+c_{\tau}\partial_{t}L+c=0,\\\ L(0)=l(0),\quad
L(T)=l(T),\end{cases}\quad\text{where }c_{\tau}=\begin{cases}c\quad\text{on
}(0,\tau),\\\ -c\quad\text{on }(\tau,T),\end{cases}$

for some $\tau\in(0,T)$, depending on the sign of $\partial_{t}L$. Let us
first consider the problem

$\begin{cases}\partial_{t}^{2}v+c_{\tau}\partial_{t}v=0,\\\ v(0)=l(0),\quad
v(T)=l(T).\end{cases}$

Observe that, since $\partial_{t}v$ satisfies a homogeneous linear
differential equation of first order, it does not change its sign, which
depends only on the values $v(0)$ and $v(T)$. Then, by the maximum principle
we can say that $v\leq V$ on $[0,T]$, with $V$ solution to

$\begin{cases}\partial_{t}^{2}V+c\partial_{t}V=0,\\\ V(0)=l(0),\quad
V(T)=l(T).\end{cases}$

Computing directly the solution $V$ we find that

$v(t)\leq V(t)=l(0)\frac{e^{-\gamma t}-e^{-\gamma T}}{1-e^{-\gamma
T}}+l(T)\frac{1-e^{-\gamma t}}{1-e^{-\gamma T}},$

for any $t\in[0,T]$, with $\gamma$ being either $c$ or $-c$, depending on
whether $l(0)\leq l(T)$ or $l(0)>l(T)$. Then, one can also find an explicit
solution $w$ of the non-homogeneous equation with zero boundary data and show
that

$w(t)\leq\frac{2e^{cT}}{c}\quad\text{for any }t\in(0,T).$

We refer to [33, Lemma 3.1.5] for all these explicit calculations. Hence, we
find that

$l(t)\leq L(t)=v(t)+w(t)\leq l(0)\frac{e^{-\gamma t}-e^{-\gamma
T}}{1-e^{-\gamma T}}+l(T)\frac{1-e^{-\gamma t}}{1-e^{-\gamma
T}}+\frac{2e^{cT}}{c},$

for any $t\in[0,T]$, with $\gamma=c$ or $\gamma=-c$. Therefore, by using the
definition of $l$ and the fact that $\int_{0}^{t}g\,\mathrm{d}s\leq 2CT$, we
deduce that

$\log\lVert\boldsymbol{\uppsi}(t)\rVert\leq\log\lVert\boldsymbol{\uppsi}(0)\rVert\frac{e^{-\gamma
t}-e^{-\gamma T}}{1-e^{-\gamma
T}}+\log\lVert\boldsymbol{\uppsi}(T)\rVert\frac{1-e^{-\gamma t}}{1-e^{-\gamma
T}}+\frac{2e^{cT}}{c}+2CT.$

Finally, by exponentiating and using
$\lVert\boldsymbol{\uppsi}(0)\rVert_{\mathbb{H}}\leq M$, we reach the thesis
with a constant $C_{1}>0$ depending only on $c$, $C$ and $T$. ∎

######  Remark 3.2.

In the proof of Proposition 3.1, for simplicity of exposition, we used the
abstract formulation
$\partial_{t}\boldsymbol{\uppsi}-\boldsymbol{D}\Delta\boldsymbol{\uppsi}=f_{\boldsymbol{\uppsi}}$,
$\boldsymbol{\uppsi}(0)=\boldsymbol{\uppsi}_{0}$ of the system (3.4)–(3.8). We
just want to comment on the fact that the Laplacian operator is not exactly
the same for all three components, since for $\varphi$ it is complemented with
homogeneous Dirichlet boundary conditions, while for $\sigma$ and $p$ with
homogeneous Neumann. However, with a slight abuse of notation, we will stick
to the already-used expression.

######  Remark 3.3.

Observe that, from [33, Lemma 3.1.5], the choice of $\lambda_{1}$ or
$\lambda_{2}$ depends on whether
$\log\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}$ is smaller or bigger
than
$\log\lVert\boldsymbol{\uppsi}(T)\rVert^{2}_{\mathbb{H}}-\int_{0}^{T}(f_{\boldsymbol{\uppsi}},\boldsymbol{\uppsi})/\lVert\boldsymbol{\uppsi}\rVert_{\mathbb{H}}^{2}\,\mathrm{d}t$,
which in general cannot be verified a priori. Therefore, we need to take into
account both versions of $\lambda(t)$. However, in both cases
$\lambda_{i}(0)=0$, $i=1,2$, so (3.2) does not give a Hölder stability
estimate also for $t=0$. Nevertheless, we still get a backward uniqueness
result as in the following Corollary.

######  Corollary 3.4.

Assume hypotheses A1–A3. Let $(\varphi_{1},\sigma_{1},p_{1})$ and
$(\varphi_{2},\sigma_{2},p_{2})$ be two solutions of (2.1)–(2.5) corresponding
to two triples of initial data
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})\in\mathcal{I}_{\text{ad}}$ for
$i=1,2$.

If $(\varphi_{1},\sigma_{1},p_{1})(T)=(\varphi_{2},\sigma_{2},p_{2})(T)$, then
$(\varphi_{1},\sigma_{1},p_{1})(t)=(\varphi_{2},\sigma_{2},p_{2})(t)$ for any
$t\in[0,T]$. In particular,
$(\varphi_{0}^{1},\sigma_{0}^{1},p_{0}^{1})=(\varphi_{0}^{2},\sigma_{0}^{2},p_{0}^{2})$
in $\mathbb{V}$.

###### Proof.

Observe that, if
$(\varphi_{1},\sigma_{1},p_{1})(T)=(\varphi_{2},\sigma_{2},p_{2})(T)$, then
estimate (3.2) gives that

$\lVert(\varphi_{1}(t),\sigma_{1}(t),p_{1}(t))-(\varphi_{2}(t),\sigma_{2}(t),p_{2}(t))\rVert_{\mathbb{H}}=0$

for any $t\in(0,T]$. The result now easily follows also for $t=0$, by
exploiting the continuity of the function
$t\to\lVert(\varphi_{1}(t),\sigma_{1}(t),p_{1}(t))-(\varphi_{2}(t),\sigma_{2}(t),p_{2}(t))\rVert_{\mathbb{H}}$,
which is guaranteed by the $\mathcal{C}^{0}([0,T];\mathbb{H})$-regularity of
the solutions. ∎

For future use, we observe that if we define (cf. (3.1) for the definition of
$M$)

$\varepsilon:=\frac{\lVert(\varphi_{1}(T),\sigma_{1}(T),p_{1}(T))-(\varphi_{2}(T),\sigma_{2}(T),p_{2}(T))\rVert_{\mathbb{H}}}{M}=\frac{\lVert\boldsymbol{\uppsi}(T)\rVert_{\mathbb{H}}}{M},$
(3.14)

then the stability estimate (3.2) becomes

$\lVert(\varphi_{1}(t),\sigma_{1}(t),p_{1}(t))-(\varphi_{2}(t),\sigma_{2}(t),p_{2}(t))\rVert_{\mathbb{H}}\leq
C_{1}M\varepsilon^{\lambda(t)},$

for any $t\in(0,T)$, where $\lambda(t)$ is either $\lambda_{1}(t)$ or
$\lambda_{2}(t)$. To obtain a stability estimate even for the reconstruction
of the initial data, we have to impose a further a priori bound. Indeed, if,
additionally to $(\varphi_{0},\sigma_{0},p_{0})\in\mathcal{I}_{\text{ad}}$, we
also assume that

$\lVert(\varphi_{0},\sigma_{0},p_{0})\rVert_{\mathbb{V}}\leq\overline{C},$

for some $\overline{C}>0$, then by Theorem 2.1 we have the following uniform
bound on the solution of (2.1)–(2.5):

$\lVert(\varphi,\sigma,p)\rVert_{H^{1}(0,T;\mathbb{H})}\leq M_{1},$ (3.15)

where $M_{1}$ depends only on the fixed parameters of the system and on
$\overline{C}$, and can be chosen as the minimal one. This further bound
enables us to prove a logarithmic-type stability result for the reconstructed
initial data.

######  Theorem 3.5.

Assume hypotheses A1–A3. Let $(\varphi_{1},\sigma_{1},p_{1})$ and
$(\varphi_{2},\sigma_{2},p_{2})$ be two solutions of (2.1)–(2.5) corresponding
to two triples of initial data
$(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})\in\mathcal{I}_{\text{ad}}$ for
$i=1,2$.

Assume that for $i=1,2$ we also have the following a priori bound:

$\lVert(\varphi_{0}^{i},\sigma_{0}^{i},p_{0}^{i})\rVert_{\mathbb{V}}\leq\overline{C},$
(3.16)

and let $M>0$, $C_{1}>0$, $M_{1}>0$ be as in (3.1), (3.2), (3.15). Let
$\varepsilon$ be defined as in (3.14) and assume that

$\varepsilon\leq\exp\left\\{-\left(\min\left\\{1,\frac{4\sqrt{3}MC_{1}^{3/2}}{9M_{1}}\right\\}\right)^{-1}\right\\}.$

Then, there exists a constant $C_{2}>0$ such that

$\lVert(\varphi_{0}^{1},\sigma_{0}^{1},p_{0}^{1})-(\varphi_{0}^{2},\sigma_{0}^{2},p_{0}^{2})\rVert^{2}_{\mathbb{H}}\leq\frac{C_{2}}{\sqrt{\lvert\log\varepsilon\rvert}},$
(3.17)

where we can quantify the constant as

$C_{2}=\frac{2M_{1}MC_{1}^{1/2}}{\beta^{1/2}}+\frac{3M_{1}^{2}}{4\beta
C_{1}},\quad\text{with }\beta=\frac{\gamma}{e^{\gamma T}-1}>0,$

and where $\gamma>0$ is the constant appearing in (3.3).

###### Proof.

For all this argument, we still use the compact notation with
$\boldsymbol{\uppsi}$ introduced in the proof of Proposition 3.1. Note that
under our regularity assumptions, by standard results in Lebesgue-Bochner
theory, the function $t\to\lVert\boldsymbol{\uppsi}(t)\rVert^{2}_{\mathbb{H}}$
is absolutely continuous. Then, we can apply the fundamental theorem of
calculus to say that for any $t\in[0,T]$

$\displaystyle\lVert\boldsymbol{\uppsi}(t)\rVert^{2}_{\mathbb{H}}$
$\displaystyle=\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}+\int_{0}^{t}\frac{\mathrm{d}}{\mathrm{d}s}\lVert\boldsymbol{\uppsi}(s)\rVert^{2}_{\mathbb{H}}\,\mathrm{d}s$
$\displaystyle=\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}+\int_{0}^{t}2(\partial_{s}\boldsymbol{\uppsi}(s),\boldsymbol{\uppsi}(s))_{\mathbb{H}}\,\mathrm{d}s.$

Then, by using Cauchy-Schwarz, the stability inequality (3.2) and the bounds
(3.1) and (3.15), we deduce that

$\displaystyle\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}$
$\displaystyle\leq\lVert\boldsymbol{\uppsi}(t)\rVert^{2}_{\mathbb{H}}+2\left(\int_{0}^{t}\lVert\partial_{s}\boldsymbol{\uppsi}(s)\rVert^{2}_{\mathbb{H}}\,\mathrm{d}s\right)^{1/2}\left(\int_{0}^{t}\lVert\boldsymbol{\uppsi}(s)\rVert^{2}_{\mathbb{H}}\,\mathrm{d}s\right)^{1/2}$
$\displaystyle\leq
C_{1}^{2}M^{2}\varepsilon^{2\lambda(t)}+2M_{1}\left(\int_{0}^{t}M\cdot
C_{1}M\varepsilon^{\lambda(s)}\,\mathrm{d}s\right)^{1/2},$

where $\lambda$ can be either $\lambda_{1}$ or $\lambda_{2}$. Now, by working
on the explicit expressions of $\lambda_{i}(t)$, $i=1,2$, one can easily see
that $\lambda_{i}(0)=0$, $\lambda_{i}(T)=1$, $\lambda_{i}$ is strictly
increasing for both $i=1,2$, $\lambda_{1}$ is concave and $\lambda_{2}$ is
convex. So one has the situation depicted in Figure 3. Then, the tangent line
at $0$ to $\lambda_{2}$, which has the expression

$\tilde{\lambda}(t)=\beta t,\quad\text{with }\beta=\frac{\gamma}{e^{\gamma
T}-1},$

is below both $\lambda_{1}$ and $\lambda_{2}$. Therefore, since
$0<\varepsilon<1$, we can infer that

$\varepsilon^{\lambda_{i}(t)}\leq\varepsilon^{\beta t}\quad\text{for both
}i=1,2.$

Hence, we can further estimate

$\displaystyle\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}$
$\displaystyle\leq C_{1}^{2}M^{2}\varepsilon^{2\beta
t}+2M_{1}MC^{1/2}\left(\int_{0}^{t}\varepsilon^{\beta
s}\,\mathrm{d}s\right)^{1/2}$ $\displaystyle=C_{1}^{2}M^{2}\varepsilon^{2\beta
t}+2M_{1}MC^{1/2}\left(\frac{\varepsilon^{\beta
t}}{\beta\log\varepsilon}-\frac{1}{\beta\log\varepsilon}\right)^{1/2}$
$\displaystyle=C_{1}^{2}M^{2}\varepsilon^{2\beta
t}+\frac{2M_{1}MC^{1/2}}{\beta^{1/2}\sqrt{\lvert\log\varepsilon\rvert}}(1-\varepsilon^{\beta
t})^{1/2}:=g(t).$

Now, note that this estimate holds for any $t\in[0,T]$, so we can make it
sharper and independent of $t$ by minimising the function $g(t)$. Then, we
compute its derivative and, recalling that $\log\varepsilon<0$, we obtain that

$g^{\prime}(t)=-C_{1}^{2}M^{2}2\beta\lvert\log\varepsilon\rvert\cdot\varepsilon^{2\beta
t}+M_{1}MC_{1}^{1/2}\beta^{1/2}\sqrt{\lvert\log\varepsilon\rvert}\frac{\varepsilon^{\beta
t}}{(1-\varepsilon^{\beta t})^{1/2}}.$

If we call $x=\varepsilon^{\beta t}$, $0<x<1$, we can find the critical points
by solving the equation

$-C_{1}^{2}M^{2}2\beta\lvert\log\varepsilon\rvert\cdot
x^{2}+M_{1}MC_{1}^{1/2}\beta^{1/2}\sqrt{\lvert\log\varepsilon\rvert}\frac{x}{(1-x)^{1/2}}=0,$

which, after some manipulations, is equivalent to

$x\sqrt{1-x}=\frac{M_{1}}{2\beta^{1/2}C_{1}^{3/2}M\sqrt{\lvert\log\varepsilon\rvert}}.$
(3.18)

By analysing the function $f(x)=x\sqrt{1-x}$, we see that its graph is like
the one in Figure 4, with maximum value at $\frac{2\sqrt{3}}{9}$. Then,
equation (3.18) has at least a solution if

$\frac{M_{1}}{2\beta^{1/2}C_{1}^{3/2}M\sqrt{\lvert\log\varepsilon\rvert}}\leq\frac{2\sqrt{3}}{9},$

which gives part of the smallness condition in the statement of the Theorem.
Moreover, one can easily see that the solution corresponding to the minimum
value of $g(t)$ is the smallest one, let us call it $\overline{x}$. Then, we
can see that the increasing branch of the function $f$ lies between the two
lines $y=x$ and $y=\frac{\sqrt{3}}{3}x$, therefore we can estimate the value
of $\overline{x}$ with the solutions corresponding to the two lines, namely

$\frac{M_{1}}{2\beta^{1/2}C_{1}^{3/2}M\sqrt{\lvert\log\varepsilon\rvert}}\leq\overline{x}\leq\frac{\sqrt{3}M_{1}}{2\beta^{1/2}C_{1}^{3/2}M\sqrt{\lvert\log\varepsilon\rvert}}.$
(3.19)

Finally, we can use this piece of information to close our estimate:

$\displaystyle\lVert\boldsymbol{\uppsi}(0)\rVert^{2}_{\mathbb{H}}$
$\displaystyle\leq\min_{t\in[0,T]}g(t)=C_{1}^{2}M^{2}\overline{x}^{2}+\frac{2M_{1}MC^{1/2}}{\beta^{1/2}\sqrt{\lvert\log\varepsilon\rvert}}(1-\overline{x})^{1/2}$
$\displaystyle\leq C_{1}^{2}M^{2}\frac{3M_{1}^{2}}{4\beta
C_{1}^{3}M^{2}\lvert\log\varepsilon\rvert}+\frac{2M_{1}MC^{1/2}}{\beta^{1/2}\sqrt{\lvert\log\varepsilon\rvert}}\underbrace{\left(1-\frac{M_{1}}{2\beta^{1/2}C_{1}^{3/2}M\sqrt{\lvert\log\varepsilon\rvert}}\right)^{1/2}}_{\leq
1}$
$\displaystyle\leq\frac{1}{\sqrt{\lvert\log\varepsilon\rvert}}\left(\frac{2M_{1}MC^{1/2}}{\beta^{1/2}}+\frac{3M_{1}^{2}}{4\beta
C_{1}\sqrt{\lvert\log\varepsilon\rvert}}\right)$
$\displaystyle\leq\frac{1}{\sqrt{\lvert\log\varepsilon\rvert}}\left(\frac{2M_{1}MC^{1/2}}{\beta^{1/2}}+\frac{3M_{1}^{2}}{4\beta
C_{1}}\right),$

if we also assume that $\frac{1}{\sqrt{\lvert\log\varepsilon\rvert}}\leq 1$,
which gives the other part of the smallness assumption. This concludes the
proof of Theorem 3.5. ∎

Figure 3: Graphs of $\lambda_{1}(t)$, $\lambda_{2}(t)$ and
$\tilde{\lambda}(t)$,
where we set $T=1$ for convenience.

Figure 4: Graph of the function $f(x)=x\sqrt{1-x}$,
along with the two linear bounds.

######  Remark 3.6.

Observe that this type of logarithmic estimate can be impractical in
applications, as the error actually gets small only if the final data are very
close to each other. Moreover, it can be easily seen that the stability
constant $C_{2}$ deteriorates exponentially with the final time $T$, making
this stability estimate even more unreliable. Nevertheless, this is expected
since we are dealing with a backward problem for a parabolic system, which is
known to be severely ill-posed as $T$ gets larger.

We now want to improve this stability estimate to one of Lipschitz type, under
the additional a priori assumption that the initial data to be reconstructed
actually lives in a finite-dimensional subspace of $\mathbb{V}$, for instance,
one of the discrete spaces used in numerical approximations. To do this, we
recall the following general result, which was proved in [7, Proposition 5].
We note that the same result was successfully used to get Lipschitz stability
estimates for other inverse ill-posed problems, like in [3, 6, 10].

######  Lemma 3.7.

Let $R$ and $r_{0}$ be positive numbers. Let $\Lambda$ be an open subset of
$\mathbb{R}^{n}$, for some $n\in\mathbb{N}$, and let $K$ be a subset of
$\Lambda$. Assume that

$K\subseteq B_{R}(0)\quad\text{and}\quad B_{r_{0}}(p)\subseteq\Lambda,\text{
for any }p\in K.$ (3.20)

Let $\mathcal{B}$ be a Banach space and $G:\Lambda\to\mathcal{B}$ an operator.
Assume that $G$ satisfies the following conditions:

  * $(i)$

$G_{\mid K}$ is injective,

  * $(ii)$

$(G_{\mid K})^{-1}:G(K)\to K$ is uniformly continuous with modulus of
continuity $\omega_{0}(\cdot)$,

  * $(iii)$

$G$ is Fréchet-differentiable with derivative $\mathrm{D}G$,

  * $(iv)$

$\mathrm{D}G:\Lambda\to\mathcal{L}(\Lambda,\mathcal{B})$ is uniformly
continuous with modulus of continuity $\omega_{1}(\cdot)$.

  * $(v)$

there exists $m_{0}>0$ such that

$\inf_{x\in K,\lvert y\rvert=1}\lVert\mathrm{D}G(x)[y]\rVert_{\mathcal{B}}\geq
m_{0}.$

Then, for any $x_{1},x_{2}\in K$

$\lvert x_{1}-x_{2}\rvert\leq C\lVert G(x_{1})-G(x_{2})\rVert_{\mathcal{B}},$
(3.21)

where
$C=\max\left\\{\frac{2R}{(\omega_{0}^{-1})_{*}(\delta_{1})},\frac{2}{m_{0}}\right\\}$,
$\delta_{1}=\frac{1}{2}\min\\{\delta_{0},r_{0}\\}$ and
$\delta_{0}=(\omega_{1}^{-1})_{*}\left(\frac{m_{0}}{2}\right)$. Here, for a
modulus of continuity $\omega$, we denote

$(\omega^{-1})_{*}(y):=\inf\\{x\in[0,+\infty)\mid\omega(x)\geq y\\}.$

Before using Lemma 3.7 in our case, we still need to study more in detail the
Fréchet-derivative of our operator $\mathcal{R}:\mathbb{V}\to\mathbb{H}$,
which associates to any initial data $(\varphi_{0},\sigma_{0},p_{0})$ the
solution $(\varphi(T),\sigma(T),p(T))$ at the final time. In Theorem 2.5 we
showed that $\mathrm{D}\mathcal{R}:\mathbb{V}\cap
L^{\infty}(\Omega)^{3}\to\mathcal{L}(\mathbb{V}\cap
L^{\infty}(\Omega)^{3},\mathbb{H})$ is such that, for any
$(h,k,w)\in\mathbb{V}\cap L^{\infty}(\Omega)^{3}$,
$\mathrm{D}\mathcal{R}(\varphi_{0},\sigma_{0},p_{0})[h,k,w]=(Y(T),Z(T),P(T))$,
which is the solution of the linearised system (2.21)–(2.25) evaluated at the
final time. Regarding this linearised inverse problem, we can prove a similar
logarithmic stability estimate, since it essentially has the same formal
structure as (3.4)–(3.8).

######  Theorem 3.8.

Assume hypotheses A1–A3. Let
$(\overline{\varphi},\overline{\sigma},\overline{p})$ be a solution of
(2.1)–(2.5) corresponding to a triple of initial data
$(\varphi_{0},\sigma_{0},p_{0})\in\mathcal{I}_{\text{ad}}$. Moreover, let
$(Y,Z,P)$ be a solution of (2.21)–(2.25) corresponding to
$(h,k,w)\in\mathbb{V}$. Assume also that

$\lVert(h,k,w)\rVert_{\mathbb{V}}\leq\overline{Q},$

for some $\overline{Q}>0$ and call $L=L(\overline{Q})>0$,
$L_{1}=L_{1}(\overline{Q})>0$ the minimal constants such that by Proposition
_2.4_

$\lVert(Y,Z,P)\rVert_{\mathcal{C}^{0}([0,T];\mathbb{H})}\leq
L,\quad\lVert(Y,Z;P)\rVert_{H^{1}(0,T;\mathbb{H})}\leq L_{1}.$ (3.22)

Then, there exists a constant $Q_{1}>0$, depending only on the parameters of
the system, such that we have the following conditional stability estimate for
any $t\in[0,T]$

$\lVert(Y(t),Z(t),P(t))\rVert_{\mathbb{H}}\leq
Q_{1}L^{1-\lambda}\lVert(Y(T),Z(T),P(T))\rVert^{\lambda},$

where $\lambda$ is either $\lambda_{1}(t)$ or $\lambda_{2}(t)$.

Moreover, if we additionally assume that

$\varepsilon:=\frac{\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}}{L}\leq\exp\left\\{-\left(\min\left\\{1,\frac{4\sqrt{3}LQ_{1}^{3/2}}{9L_{1}}\right\\}\right)^{-1}\right\\},$

then there exists a constant $Q_{2}>0$ such that

$\lVert(h,k,w)\rVert^{2}_{\mathbb{H}}\leq\frac{Q_{2}}{\sqrt{\lvert\log\varepsilon\rvert}},$
(3.23)

where we can quantify the constant as

$Q_{2}=\frac{2L_{1}LC_{1}^{1/2}}{\beta^{1/2}}+\frac{3L_{1}^{2}}{4\beta
Q_{1}}.$ (3.24)

######  Remark 3.9.

Since the system is linear, it is clear that (3.23) gives a stability estimate
for the linearised inverse problem of reconstructing the initial data
$(h,k,w)$ from the final measurements $(Y(T),Z(T),P(T))$; indeed, the
difference of two solutions still satisfies the same linear system.

###### Proof.

Observe that we can rewrite the linearised system (2.21)–(2.25) in the
abstract form:

$\partial_{t}\boldsymbol{Y}-\boldsymbol{D}\Delta\boldsymbol{Y}=f_{\boldsymbol{Y}},\quad\boldsymbol{Y}(0)=\boldsymbol{h},$

where $\boldsymbol{Y}=(Y,Z,P)^{\top}$, $\boldsymbol{h}=(h,k,w)^{\top}$ and

$f_{\boldsymbol{Y}}=\begin{pmatrix}-F^{\prime\prime}(\overline{\varphi})Y+m(\overline{\sigma})h^{\prime\prime}(\overline{\varphi})Y+m^{\prime}(\overline{\sigma})h^{\prime}(\overline{\varphi})Z\\\
-\gamma_{h}Z-\gamma_{ch}\overline{\sigma}Y-\gamma_{ch}\overline{\varphi}Z+S_{ch}Y\\\
-\gamma_{p}P+\alpha_{ch}Y\end{pmatrix}.$

We recall that the Laplacian operator actually depends on the different
boundary conditions. By using the fact that
$\overline{\varphi},\overline{\sigma}\in L^{\infty}(Q_{T})$ by Theorem 2.1, it
is easy to see that

$\lVert f_{\boldsymbol{Y}}\rVert_{\mathbb{H}}\leq
C\lVert\boldsymbol{Y}\rVert_{\mathbb{H}}.$

Therefore, we are exactly in the same situation to apply [33, Theorem 3.1.3]
and then follow the same reasoning of the proof of Theorem 3.5. This concludes
the proof. ∎

We are now in the position to apply Lemma 3.7 and establish a Lipschitz
stability estimate, when starting from a finite-dimensional subspace.

######  Theorem 3.10.

Assume hypotheses A1–A3. Let $\Lambda$ be a finite-dimensional subspace of
$\mathbb{V}$ and $K\subseteq\Lambda$ such that

$K=\\{\boldsymbol{u}\in\Lambda\mid\lVert\boldsymbol{u}\rVert_{\mathbb{V}}\leq\overline{C}\\},$

for some constant $\overline{C}>0$ to be decided a priori. Assume further that

$(\varphi_{0},\sigma_{0},p_{0})\in K\cap\mathcal{I}_{\text{ad}}.$

Then, there exists a constant $C_{s}>0$ such that for any choice of
$(\varphi_{0}^{1},\sigma_{0}^{1},p_{0}^{1})$ and
$(\varphi_{0}^{2},\sigma_{0}^{2},p_{0}^{2})$ in $K\cap\mathcal{I}_{\text{ad}}$
the following stability estimate holds

$\lVert(\varphi_{0}^{1},\sigma_{0}^{1},p_{0}^{1})-(\varphi_{0}^{2},\sigma_{0}^{2},p_{0}^{2})\rVert_{\mathbb{H}}\leq
C_{s}\lVert(\varphi_{1}(T),\sigma_{1}(T),p_{1}(T))-(\varphi_{2}(T),\sigma_{2}(T),p_{2}(T))\rVert_{\mathbb{H}},$
(3.25)

where the constant $C_{s}$ can be quantified as

$C_{s}=\max\left\\{\frac{2\overline{C}}{M}e^{\frac{16C_{0}^{2}C_{2}}{m_{0}}},\frac{2}{m_{0}}\right\\},\text{
with }m_{0}=\frac{L}{C_{\Lambda}}e^{-Q^{2}_{2}}\text{ and
}C_{\Lambda}=\sup_{\boldsymbol{h}\in\Lambda\setminus\\{0\\}}\frac{\lVert\boldsymbol{h}\rVert_{\mathbb{V}}}{\lVert\boldsymbol{h}\rVert_{\mathbb{H}}},$

where $C_{0}$ is the Lipschitz constant of the Fréchet-derivative, given by
(2.28), and $M$, $L$, and $Q_{2}$ are defined in (3.1), (3.22), (3.24).

###### Proof.

We want to apply Lemma 3.7, so we have to verify its hypotheses. Indeed, we
set $X=\mathbb{V}$ and $\Lambda$ as above. Being a finite-dimensional subspace
of $\mathbb{V}$ means that $\Lambda=\mathbb{R}^{n}$ in the language of Lemma
3.7, therefore we can immediately see that $K$ satisfies hypothesis (3.20). We
consider the operator $\mathcal{R}:\Lambda\subseteq\mathbb{V}\to\mathbb{H}$
defined as
$\mathcal{R}((\varphi_{0},\sigma_{0},p_{0}))=(\varphi(T),\sigma(T),p(T))$.
Hypothesis $(i)$ on the injectivity of $\mathcal{R}_{\mid K}$ follows directly
from Corollary 3.4. Also hypothesis $(ii)$ follows immediately from Theorem
3.5, with modulus of continuity
$\omega_{0}(r)=C_{2}/\sqrt{\lvert\log(r/M)\rvert}$. Then, hypotheses $(iii)$
and $(iv)$ on Fréchet-differentiability follow from Theorem 2.5, with modulus
of continuity $\omega_{1}(r)=C_{0}r$. Finally, we have to check condition
$(v)$. To do this, notice that for any $(\varphi_{0},\sigma_{0},p_{0})\in K$
and any $(h,k,w)\in\Lambda$, we have that
$\mathrm{D}\mathcal{R}(\varphi_{0},\sigma_{0},p_{0})[h,k,w]=(Y(T),Z(T),P(T))$,
where $(Y,Z,T)$ is the solution to (2.21)–(2.25). Now, by Theorem 3.8, we can
use (3.23) to say that, if $\lVert(h,k,w)\rVert_{\mathbb{V}}\leq\overline{Q}$,
it holds that

$\lVert(h,k,w)\rVert^{2}_{\mathbb{H}}\leq\frac{Q_{2}}{\sqrt{\left\lvert\log\frac{\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}}{L}\right\rvert}},$

which implies that

$\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}\geq
Le^{-\frac{Q_{2}^{2}}{\lVert(h,k,w)\rVert^{4}_{\mathbb{H}}}}.$

Hence, we can verify condition $(v)$ as

$\displaystyle\inf_{\begin{subarray}{c}(\varphi_{0},\sigma_{0},p_{0})\in K\\\
(h,k,w)\in\Lambda\setminus\\{0\\}\end{subarray}}\frac{\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}}{\lVert(h,k,w)\rVert_{\mathbb{V}}}\geq\inf_{\begin{subarray}{c}(\varphi_{0},\sigma_{0},p_{0})\in
K\\\
(h,k,w)\in\Lambda\setminus\\{0\\}\end{subarray}}\frac{\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}}{C_{\Lambda}\lVert(h,k,w)\rVert_{\mathbb{H}}}$
$\displaystyle\quad=\inf_{\begin{subarray}{c}(\varphi_{0},\sigma_{0},p_{0})\in
K\\\
\lVert(h,k,w)\rVert_{\mathbb{H}}=1\end{subarray}}\frac{\lVert(Y(T),Z(T),P(T))\rVert_{\mathbb{H}}}{C_{\Lambda}\lVert(h,k,w)\rVert_{\mathbb{H}}}\geq\frac{L}{C_{\Lambda}}e^{-Q_{2}^{2}},$

where we used the linearity of $\mathrm{D}\mathcal{R}$ and, since $\Lambda$ is
finite-dimensional and thus all norms on it are equivalent, we estimated
$\lVert(h,k,w)\rVert_{\mathbb{V}}\leq
C_{\Lambda}\lVert(h,k,w)\rVert_{\mathbb{H}}$.

Then, we can apply Lemma 3.7 and conclude the proof of Theorem 3.10. ∎

######  Remark 3.11.

We would like to point out that the constant $C_{s}$ blows up exponentially as
the dimension of $\Lambda$ goes to infinity. Indeed, as one can see from its
expression, $C_{s}$ has a direct exponential proportionality to the constant
$C_{\Lambda}$, which is finite only if $\Lambda$ has a finite dimension and
blows up as it becomes larger. Additionally, one can also notice that the
dependence of $C_{s}$ on the final time is even worse, namely $\log(C_{s})$
depends exponentially on $T$. This means that, even if we have some kind of
Lipschitz stability, one has to be very careful when designing numerical
algorithms to approximate the solution.

##  4 Conclusions

In this paper we consider a phase field model for prostate cancer growth
previously introduced and studied in [13, 16]. While improving previous well-
posedness results, we rigorously analyse an inverse problem aiming to
reconstruct an earlier tumour state starting from a single measurement at a
later time. We first derive a conditional logarithmic stability estimate and
we subsequently prove that the solution to the inverse problem is unique.
Furthermore, by assuming to know a priori that the solution lies in a finite-
dimensional subspace, we are also able to obtain an optimal Lipschitz
stability estimate, as in Theorem 3.10. A stability estimate of this kind
generally opens the possibility of using iterative reconstruction algorithms,
such as the Landweber scheme, to approximate numerically the solution to the
inverse problem. Indeed, under the validity of a Lipschitz stability estimate
and the Fréchet-differentiability of the forward map (or more generally of the
so-called tangential cone condition), it is shown that such iterative methods
are locally convergent [17, 27, 36, 30]. Therefore, they should provide
faithful reconstructions of the earlier tumour state (i.e., the initial
conditions of the forward problem), at least if the final time $T$ is small.
On the other hand, due to the severe ill-posedness of the non-linear inverse
problem, as shown in Theorem 3.10, the stability constants become doubly
exponentially large with $T$. This means that numerical analysis of this
problem for large final times $T$ is extremely challenging. However, in our
specific cancer growth application, it is important to be able to efficiently
treat even the case of long time horizons, to better understand how the tumour
has developed from a small mass to the situation observed in the medical
images at diagnosis. This information on tumour origin and early dynamics
holds important value to estimate prognosis and survival, anticipate potential
complications derived from disease progression, and contribute to the planning
of effective treatments for each individual patient [41, 12, 19, 35, 8, 4].
For this reason, the objective of a future part of our project is to carry out
an in-depth simulation study to validate our analytical findings and propose
efficient reconstruction methods.

Acknowledgements. C. Cavaterra, M. Fornoni and E. Rocca have been partially
supported by the MIUR-PRIN Grant 2020F3NCPX “Mathematics for industry 4.0
(Math4I4)”. C. Cavaterra has been partially supported by the MIUR-PRIN Grant
2022 “Partial differential equations and related geometric-functional
inequalities”. E. Rocca also acknowledges the support of Next Generation EU
Project No.P2022Z7ZAJ (A unitary mathematical framework for modelling muscular
dystrophies). The research of E. Beretta has been supported by the Project
AD364 - Fund of NYU Abu Dhabi. The research of C. Cavaterra is part of the
activities of “Dipartimento di Eccellenza 2023-2027” of Università degli Studi
di Milano. C. Cavaterra, M. Fornoni and E. Rocca are members of GNAMPA (Gruppo
Nazionale per l’Analisi Matematica, la Probabilità e le loro Applicazioni) of
INdAM (Istituto Nazionale di Alta Matematica). G. Lorenzo acknowledges the
support of a fellowship from ”la Caixa” Foundation (ID 100010434). The
fellowship code is LCF/BQ/PI23/11970033.

## References

  * [1] S. Agmon and L. Nirenberg, _Lower bounds and uniqueness theorems for solutions of differential equations in a Hilbert space_ , Comm. Pure Appl. Math. 20 (1967), 207–229.
  * [2] V. Akcelik, G. Biros, A. Draganescu, O. Ghattas, J. Hill, and B. van Bloemen Waanders, _Inversion of airborne contaminants in a regional model_ , Computational Science–ICCS 2006: 6th International Conference, Reading, UK, May 28-31, 2006. Proceedings, Part III 6, Springer, 2006, pp. 481–488.
  * [3] G. Alessandrini and S. Vessella, _Lipschitz stability for the inverse conductivity problem_ , Adv. in Appl. Math. 35 (2005), 207–241.
  * [4] A. Ali, T. Elumalai, B. Venkatesulu, L. Hekman, H. Mistry, et al., _Tale of two zones: investigating the clinical outcomes and research gaps in peripheral and transition zone prostate cancer through a systematic review and meta-analysis_ , BMJ Oncology 3 (2024), e000193.
  * [5] P. Arratia, M. Courdurier, E. Cueva, A. Osses, and B. Palacios, _Lipschitz stability for backward heat equation with application to fluorescence microscopy_ , SIAM J. Math. Anal. 53 (2021), 5948–5978.
  * [6] A. Aspri, E. Beretta, E. Francini, and S. Vessella, _Lipschitz stable determination of polyhedral conductivity inclusions from local boundary measurements_ , SIAM J. Math. Anal. 54 (2022), 5182–5222.
  * [7] V. Bacchelli and S. Vessella, _Lipschitz stability for a stationary $2{D}$ inverse problem with unknown polygonal boundary_, Inverse Problems 22 (2006), 1627–1658.
  * [8] T. Barrett, A. Rajesh, A. Rosenkrantz, P. Choyke, and B. Turkbey, _PI-RADS version 2.1: one small step for prostate MRI_ , Clinical radiology 74 (2019), 841–852.
  * [9] E. Beretta and E. Francini, _Lipschitz stability for the electrical impedance tomography problem: the complex case_ , Comm. Partial Differential Equations 36 (2011), 1723–1749.
  * [10] E. Beretta, E. Francini, and S. Vessella, _Determination of a linear crack in an elastic body from boundary measurements—Lipschitz stability_ , SIAM J. Math. Anal. 40 (2008), 984–1002.
  * [11] C. Cavaterra, E. Rocca, and H. Wu, _Long-time dynamics and optimal control of a diffuse interface model for tumor growth_ , Appl. Math. Optim. 83 (2021), 739–787.
  * [12] A. Chaudhuri, G. Pash, D. A. Hormuth, G. Lorenzo, M. Kapteyn, C. Wu, E. A. Lima, T. E. Yankeelov, and K. Willcox, _Predictive digital twin for optimizing patient-specific radiotherapy regimens under uncertainty in high-grade gliomas_ , Frontiers in Artificial Intelligence 6 (2023), 1222612\.
  * [13] P. Colli, H. Gomez, G. Lorenzo, G. Marinoschi, A. Reali, and E. Rocca, _Mathematical analysis and simulation study of a phase-field model of prostate cancer growth with chemotherapy and antiangiogenic therapy effects_ , Math. Models Methods Appl. Sci. 30 (2020), 1253–1295.
  * [14] P. Colli, G. Gilardi, and D. Hilhorst, _On a Cahn-Hilliard type phase field system related to tumor growth_ , Discrete Contin. Dyn. Syst. 35 (2015), 2423–2442.
  * [15] P. Colli, G. Gilardi, E. Rocca, and J. Sprekels, _Optimal distributed control of a diffuse interface model of tumor growth_ , Nonlinearity 30 (2017), 2518–2546.
  * [16] P. Colli, H. Gomez, G. Lorenzo, G. Marinoschi, A. Reali, and E. Rocca, _Optimal control of cytotoxic and antiangiogenic therapies on prostate cancer growth_ , Math. Models Methods Appl. Sci. 31 (2021), 1419–1468.
  * [17] M. V. de Hoop, L. Qiu, and O. Scherzer, _Local analysis of inverse problems: Hölder stability and iterative reconstruction_ , Inverse Problems 28 (2012), 045001, 16 pp.
  * [18] H. W. Engl and W. Rundell (eds.), _Inverse problems in diffusion processes_ , SIAM, Philadelphia, PA; GAMM, Regensburg, 1995.
  * [19] I. Y. Eyüpoglu, M. Buchfelder, and N. E. Savaskan, _Surgical resection of malignant gliomas—role in optimizing patient outcome_ , Nature Reviews Neurology 9 (2013), 141–151.
  * [20] M. Fornoni, _Maximal regularity and optimal control for a non-local Cahn-Hilliard tumour growth model_ , preprint arXiv:2402.04204 (2024).
  * [21] S. Frigeri, M. Grasselli, and E. Rocca, _On a diffuse interface model of tumour growth_ , European J. Appl. Math. 26 (2015), 215–243.
  * [22] S. Frigeri, K. F. Lam and A. Signori, _Strong well-posedness and inverse identification problem of a non-local phase field tumour model with degenerate mobilities_ , European J. Appl. Math. 33 (2022), 267–308.
  * [23] H. Garcke, K. F. Lam, and E. Rocca, _Optimal control of treatment time in a diffuse interface model of tumor growth_ , Appl. Math. Optim. 78 (2018), 495–544.
  * [24] H. Garcke, K. F. Lam, E. Sitka, and V. Styles, _A Cahn-Hilliard-Darcy model for tumour growth with chemotaxis and active transport_ , Math. Models Methods Appl. Sci. 26 (2016), 1095–1148.
  * [25] F. Giganti, A. Stabile, V. Stavrinides, E. Osinibi, A. Retter, et al., _Natural history of prostate cancer on active surveillance: stratification by MRI using the precise recommendations in a UK cohort_ , European Radiology 31 (2021), 1644–1655.
  * [26] R. K. Gupta, A. M. Roy, A. Gupta, K. Takabe, A. Dhakal, M. Opyrchal, P. Kalinski, and S. Gandhi, _Systemic therapy de-escalation in early-stage triple-negative breast cancer: dawn of a new era?_ , Cancers 14 (2022), 1856.
  * [27] M. Hanke, A. Neubauer, and O. Scherzer, _A convergence analysis of the Landweber iteration for nonlinear ill-posed problems_ , Numer. Math. 72 (1995), 21–37.
  * [28] D. N. Hào and N. V. Duc, _Stability results for backward parabolic equations with time-dependent coefficients_ , Inverse Problems 27 (2011), 025003, 20.
  * [29] D. A. Hormuth, K. A. Al Feghali, A. M. Elliott, T. E. Yankeelov, and C. Chung, _Image-based personalization of computational models for predicting response of high-grade glioma to chemoradiation_ , Sci. Rep. 11 (2021), 8520.
  * [30] S. Hubmer and R. Ramlau, _Convergence analysis of a two-point gradient method for nonlinear ill-posed problems_ , Inverse Problems 33 (2017), 095004.
  * [31] O. Y. Imanuvilov and M. Yamamoto, _Lipschitz stability in inverse parabolic problems by the Carleman estimate_ , Inverse Problems 14 (1998), 1229–1245.
  * [32] O. Y. Imanuvilov and M. Yamamoto, _Conditional stability in a backward parabolic system_ , Appl. Anal. 93 (2014), 2174–2198.
  * [33] V. Isakov, _Inverse problems for partial differential equations_ , third ed., Applied Mathematical Sciences, vol. 127, Springer, Cham, 2017.
  * [34] R. Jaroudi, G. Baravdish, B. T. Johansson, and F. Åström, _Numerical reconstruction of brain tumours_ , Inverse Probl. Sci. Eng. 27 (2019), 278–298.
  * [35] C. Jungk, R. Warta, A. Mock, S. Friauf, B. Hug, D. Capper, A. Abdollahi, J. Debus, M. Bendszus, A. von Deimling, et al., _Location-dependent patient outcome and recurrence patterns in idh1-wildtype glioblastoma_ , Cancers 11 (2019), 122.
  * [36] B. Kaltenbacher, A. Neubauer, and O. Scherzer, _Iterative regularization methods for nonlinear ill-posed problems_ , Walter de Gruyter, Berlin, 2008.
  * [37] A. S. Kazerouni, M. Gadde, A. Gardner, D. A. Hormuth, A. M. Jarrett, et al., _Integrating quantitative assays with biologically based mathematical modeling for predictive oncology_ , Iscience 23 (2020), 101807.
  * [38] I. Kukavica, _Log-log convexity and backward uniqueness_ , Proc. Amer. Math. Soc. 135 (2007), 2415–2421.
  * [39] M. M. Lavrent’ev, V. G. Romanov, and S. P. Shishat$\cdot$skiĭ, _Ill-posed problems of mathematical physics and analysis_ , Translations of Mathematical Monographs, vol. 64, American Mathematical Society, Providence, RI, 1986, Translated from Russian by J. R. Schulenberger.
  * [40] G. Lorenzo, J. S. Heiselman, M. A. Liss, M. I. Miga, H. Gomez, T. E. Yankeelov, A. Reali, and T. J. Hughes, _A pilot study on patient-specific computational forecasting of prostate cancer growth during active surveillance using an imaging-informed biomechanistic model_ , Cancer Research Communications 4 (2024), 617–633.
  * [41] G. Lorenzo, D. A. Hormuth II, A. M. Jarrett, E. A. Lima, S. Subramanian, G. Biros, J. T. Oden, T. J. Hughes, and T. E. Yankeelov, _Quantitative in vivo imaging to enable tumour forecasting and treatment optimization_ , Cancer, Complexity, Computation, Springer, 2022, pp. 55–97.
  * [42] G. Lorenzo, M. A. Scott, K. Tew, T. J. Hughes, Y. J. Zhang, L. Liu, G. Vilanova, and H. Gomez, _Tissue-scale, personalized modeling and simulation of prostate cancer growth_ , Proc. Natl. Acad. Sci. U.S.A. 113 (2016), E7663–E7671.
  * [43] J. S. Lowengrub, H. B. Frieboes, F. Jin, Y.-L. Chuang, X. Li, P. Macklin, S. M. Wise, and V. Cristini, _Nonlinear modelling of cancer: bridging the gap between cells and tumours_ , Nonlinearity 23 (2010), R1–R91.
  * [44] A. Miranville, E. Rocca, and G. Schimperna, _On the long time behavior of a tumor growth model_ , J. Differential Equations 267 (2019), 2616–2642.
  * [45] N. Mottet, R. C. van den Bergh, E. Briers, T. Van den Broeck, M. G. Cumberbatch, et al., _EAU-EANM-ESTRO-ESUR-SIOG guidelines on prostate cancer—2020 update. part 1: screening, diagnosis, and local treatment with curative intent_ , European Urology 79 (2021), 243–262.
  * [46] D. E. Neal, C. Metcalfe, J. L. Donovan, J. A. Lane, M. Davis, et al., _Ten-year mortality, disease progression, and treatment-related side effects in men with localised prostate cancer from the protect randomised controlled trial according to treatment received_ , European Urology 77 (2020), 320–330.
  * [47] L. E. Payne, _Improperly posed problems in partial differential equations_ , Regional Conference Series in Applied Mathematics, vol. No. 22, Society for Industrial and Applied Mathematics, Philadelphia, PA, 1975.
  * [48] S. Subramanian, K. Scheufele, M. Mehl, and G. Biros, _Where did the tumor start? An inverse solver with sparse localization for tumor growth models_ , Inverse problems 36 (2020), 045006.
  * [49] A. Tudorica, K. Y. Oh, S. Y. Chui, N. Roy, M. L. Troxell, et al., _Early prediction and evaluation of breast cancer response to neoadjuvant chemotherapy using quantitative dce-mri_ , Translational oncology 9 (2016), 8–17.
  * [50] S. Vessella, _Unique continuation properties and quantitative estimates of unique continuation for parabolic equations_ , Handbook of differential equations: evolutionary equations. Vol. V, Handb. Differ. Equ., Elsevier/North-Holland, Amsterdam, 2009, pp. 421–500.
  * [51] S. Wise, J. Lowengrub, H. Frieboes, and V. Cristini, _Three-dimensional multispecies nonlinear tumor growth-I: Model and numerical method_ , J. Theor. Biol. 253 (2008), 524–543.
  * [52] C. Wu, A. M. Jarrett, Z. Zhou, N. Elshafeey, B. E. Adrada, et al., _MRI-based digital models forecast patient-specific treatment responses to neoadjuvant chemotherapy in triple-negative breast cancer_ , Cancer Research 82 (2022), 3394–3404.
  * [53] J. Xu, G. Vilanova, and H. Gomez, _A mathematical model coupling tumor growth and angiogenesis_ , PloS One 11 (2016), e0149422.
  * [54] M. Yamamoto, _Carleman estimates for parabolic equations and applications_ , Inverse Problems 25 (2009), 123013, 75.
  * [55] A. Yin, D. J. A. Moes, J. G. van Hasselt, J. J. Swen, and H.-J. Guchelaar, _A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors_ , CPT: pharmacometrics & systems pharmacology 8 (2019), 720–737.
  * [56] M. Ziu, B. Y. Kim, W. Jiang, T. Ryken, and J. J. Olson, _The role of radiation therapy in treatment of adults with newly diagnosed glioblastoma multiforme: a systematic review and evidence-based clinical practice guideline update_ , Journal of Neuro-oncology 150 (2020), 215–267.