1. Introduction
In this work, we extend our previously proposed model from [Reference Dondl, Poh and Zeinhofer13] for bone regeneration in the presence of a bioresorbable porous scaffold to include a scaffold architecture optimization problem that allows geometry- and patient-dependent optimal scaffold designs. We prove the existence of optimal scaffold density distributions under assumptions on the model’s geometry and boundary conditions that are realistic for applications, i.e., non-smooth domains and mixed Dirichlet-Neumann boundary conditions. Furthermore, we present numerical simulations that indicate how to mitigate the negative impact of stress shielding on the bone regeneration process that appears in vivo due to the external fixation of the scaffold at the defect site, as documented in, e.g., [Reference Arabnejad, Johnston, Tanzer and Pasini2, Reference Behrens, Wirth, Windhagen, Nolte, Meyer-Lindenberg and Bouguecha3, Reference Huiskes, Weinans and Van Rietbergen23, Reference Sumner and Galante36] in the context of total hip arthroplasty or [Reference Terjesen, Nordby and Arnulf37, Reference Viateau, Guillemin, Bousson, Oudina, Hannouche, Sedel, Logeart‐Avramoglou and Petite38] for femoral defects.
Our model was previously proposed in [Reference Poh, Valainis, Bhattacharya, van Griensven and Dondl32] and analysed in [Reference Dondl, Poh and Zeinhofer13]. Its essential processes are an interplay between the mechanical and biological environment, which we model by a coupled system of Partial differential equation (PDE) and ordinary differential equation (ODE). The mechanical environment is represented by a linear elastic equation and the biological environment through diffusion equations as well as logistic ODEs, modelling signalling molecules and cells/bone, respectively. Material properties are incorporated using homogenized quantities, not resolving any scaffold microstructure. This makes the model efficient in computations and thus allows us to solve the scaffold architecture optimization problem numerically with manageable computational cost.
The article is organized as follows. In following section, we provide a brief introduction to tissue engineering for the treatment of severe bone defects, present our computational model and the corresponding scaffold optimization problem. Then, we discuss its weak formulation in Section 2, followed by the presentation of the main analytical results in Section 3 together with an outline of their proofs. We proceed with the discussion of the numerical simulations in Section 4. The detailed proofs of the main analytical results are provided in Appendix A.
1.1. Scaffold mediated bone growth
The treatment of critical-sized bone defects (
$\gt$
25 mm) and restoration of skeletal functions is challenging using current treatment options, see [Reference Nauth, Schemitsch, Norris, Nollin and Watson27]. Non-unions of the defect, i.e., when no bridging is achieved after
$\gt$
9 months and healing stagnates for 3 months, pose severe problems for patients, as discussed in [Reference Calori6], and have a relevant financial dimension. For the UK, [Reference Stewart35] estimates the healthcare costs to
$\unicode{x00A3} 320$
million annually. Further, healing prospects can be worsened by comorbidities such as diabetes, which is associated with compromised bone regeneration capacity, see [Reference Marin, Luyten, Van der Schueren, Kerckhofs and Vandamme26].
Recent research illustrates the potential of porous, bio-resorbable support structures, as temporary support. These structures are called scaffolds and are implanted in the defect site to provide stability, allow vascularization and guide new bone formation. Promising results were recently obtained in vivo and in clinical studies. We refer to [Reference Cipitria, Lange, Schell, Wagermaier, Reichert, Hutmacher, Fratzl and Duda8, Reference Paris, Götz, Hettrich, Bidan, Dunlop, Razi, Zizak, Hutmacher, Fratzl, Duda, Wagermaier and Cipitria28, Reference Petersen, Princ, Korus, Ellinghaus, Leemhuis, Herrera, Klaumünzer, Schreivogel, Woloszyk, Schmidt-Bleek, Geissler, Heschel and Duda29, Reference Pobloth, Checa, Razi, Petersen, Weaver, Schmidt-Bleek, Windolf, Tatai, Roth, Schaser, Duda and Schwabe31]. These works also indicate that material choice and scaffold design are critical variables for a successful healing outcome, yet they are at present not fully understood.
Several objectives need to be considered for the scaffold design, including (a) the pore size, porosity and shape of the microstructure, which have a significant influence on vascularization, cell proliferation and cell differentiation; (b) the mechanical properties of the scaffold need to guarantee a proper strain distribution within the defect site in order for bone to grow; (c) information concerning comorbidities of the patients that imply reduced bone growth or bone density as caused by diabetes needs to be incorporated in the scaffold design. Hence, designing scaffolds that are optimized for the patient and defect site at hand is of fundamental importance. With the possibilities of additive manufacturing, personalized scaffold designs are within reach.
So far, the question of optimal scaffold design has been dominated by trial-and-error approaches. However, this workflow is expensive and prohibits patient-specific designs. Topology optimization techniques have successfully been exploited for optimal design questions yielding scaffolds that meet elastic optimality with given porosity or fluid permeability, as shown in [Reference Challis, Guest, Grotowski and Roberts7, Reference Coelho, Hollister, Flanagan and Fernandes9, Reference Dias, Guedes, Flanagan, Hollister and Fernandes10, Reference Dondl, Poh, Rumpf and Simon12, Reference Guest and Prévost19, Reference Kang, Lin and Hollister24, Reference Lin, Kikuchi and Hollister25, Reference Wang, Xu, Zhou, Xu, Leary, Choong, Qian, Brandt and Xie39]. However, these methods usually lack the ability to resolve the time dependence of the elastic moduli in scaffold-mediated bone growth, which is integral to incorporate for an optimal scaffold design.
Based on the previous studies [Reference Dondl, Poh and Zeinhofer13, Reference Poh, Valainis, Bhattacharya, van Griensven and Dondl32], we propose a PDE-constrained optimization problem for the optimal porosity (or, equivalently, density) distribution of a scaffold. The advantages of optimized scaffolds include (a), the mitigation of stress shielding effects, which are caused by external fixation of the defect site and result in areas of low stress within the scaffold. This leads to poor bone regeneration, as mechanical stimulus is indispensable for bone growth. (b) Optimized scaffolds can be designed patient-dependently by altering the model’s parameters. It is straightforward to include a reduced bone regeneration capability due to, e.g., diabetes into the model.
Note that the model proposed in [Reference Dondl, Poh and Zeinhofer13] does not resolve the porous micro-structure of the scaffold design, but uses coarse-grained values instead, i.e., values averaged over a volume representative of the scaffold microstructure (a representative volume element, RVE), and thus optimizes only these macroscopic quantities. In a scaffold based on a unit cell design, the scaffold volume fraction (or equivalently, the porosity) changes on a larger length-scale than the unit cell design. Likewise, the other quantities of the model can be viewed as locally averaged values. We stress that this homogenization viewpoint of the model is essential for the feasibility of the optimization problem and we refer to [Reference Dondl, Poh and Zeinhofer13] for more information concerning this viewpoint. The resulting optimal scaffold porosity distributions can easily be introduced in a 3d printable, periodic microstructure based, scaffold design: one simply considers a one-parameter family of microstructures, where the thickness of struts or surfaces in the microstructure is given by the parameter. The parameter can then be chosen non-constant over the scaffold domain such that at each periodic unit cell the correct scaffold density is recovered as an average.
The main contributions of this article are as follows. We analyze the optimal scaffold design problem mathematically and prove the existence of optimal scaffold designs under realistic assumptions on the geometry of the computational domain, the boundary conditions and functional relationships in Section 3. The existence of such scaffolds is a necessary prerequisite for a successful numerical treatment of the optimization problem. In Section 4 we then investigate numerically how optimized scaffold designs mitigate the negative impacts of stress shielding. Our findings suggest that our three-dimensional optimization routine is able to successfully mitigate adverse stress shielding effects that may arise from external scaffold fixation.
1.2. The system of equations
The underlying paradigm of the model is that an interplay of the biological and the mechanical environment is responsible for bone growth where the mechanical environment is described through displacements and strains and the biological environment through bio-active molecules (signalling molecules) and different cell types. The model is a coupled system of evolution equations composed of a linear elastic equilibrium equation for every point in time, diffusion equations for the bio-active molecules and ordinary differential equations for the concentration of osteoblasts and the volume fraction of bone. The model considered here is a concrete instance of the more general model proposed in [Reference Dondl, Poh and Zeinhofer13].
By
$\Omega \subset \mathbb{R}^3$
we denote the computational domain, i.e., the bone defect site. We let
$I = [0,T]$
be a finite time interval. On
$\Omega$
we keep track of the local scaffold volume fraction, which we call
$\rho (x)$
, where
$x\in \Omega$
. Therefore, the porosity
$\theta$
is
$\theta (x) = 1 - \rho (x)$
, but we use only
$\rho$
. We do not resolve a time dependency for
$\rho$
. This is due to the experimental findings in [Reference Pitt, Chasalow, Hibionada, Klimas and Schindler30], which have shown that, in the time-frame of two years, PCL degrades largely due to bulk erosion. Only the molecular mass decreases, which we model by the exponential decay
$\sigma (t) = e^{-k_1 t}$
. Thus, product
$\rho (x)\cdot \sigma (t)$
quantifies the mechanical properties of PCL over time and space. We denote the bone volume fraction averaged over an RVE by
$b(t,x)$
and the variables
$b, \sigma$
and
$\rho$
then determine the material properties of the bone-scaffold composite. We work in the linear elastic regime and use the notation
$\mathbb{C}(\rho ,\sigma ,b)$
for the elastic tensor capturing the material properties.
Besides isotropy and ellipticity, we assume little for the tensor
$\mathbb{C}(\rho ,\sigma ,b)$
. Choosing a microstructure allows us to explicitly specify
$\mathbb{C}(\rho , \sigma , b)$
in numerical applications. By
$u(t,x)$
we denote the displacement field satisfying mechanical equilibrium equations, see (1). We denote the strain by
$\varepsilon (u)$
.
We represent the biological environment through two bio-active molecules
$a_1(t,x)$
and
$a_2(t,x)$
, which should be viewed as endogenous angiogenic and osteoinductive factors. These molecules diffuse depending on the scaffold density
$\rho$
which we model by
$D(\rho)$
in the equation (2). We leave this as an abstract functional relationship for the same reasons as discussed in the context of the elastic tensor. Further, we assume exponential decay and production of bio-active molecules in the presence of strain and a local density of osteoblast cells, which are denoted by
$c(t,x)$
. Essential for the production of growth factors is mechanical stimulus
$S(\varepsilon (u))$
. The quantity
$S(\varepsilon (u))$
is derived from strain; admissible choices include the strain magnitude or octahedral shear strains. We allow flexibility in the choice of
$S$
; see also the discussion in Section 2. The strain-dependent reaction term is motivated by Wolff’s law for bone remodelling, see [Reference Wolff40]. Strain as a driving force for bone regeneration is also supported by more recent work, for example in [Reference Ruff, Holt and Trinkaus33]. Note that the concentration of bio-active molecules is normalized to unity in healthy tissue. This dictates the decay and production rates in concrete simulations.
The equation (3) governs the production of osteoinductive cells (here: osteoblasts) and is modelled by logistic growth with a driving factor depending on both bio-active molecules
$a_1$
and
$a_2$
and a proliferation term
$(1+k_7c)$
. The factor
$1-c(1-\rho)^{-1}$
encodes the osteoblast carrying capacity and implies that the osteoblast density is bounded by the amount of available space, i.e., the space not filled by the scaffold. We neglect diffusion in this equation as osteoblasts diffuse on a significantly lower level than the bio-active molecules. Modelling only one cell type in our present model is a simplification and an extension of the model is easily feasible. The equation for bone growth (4) is similar to the one for osteoblast concentration. We remark that osteoblast and bone do not compete for space; this encodes the assumption that, e.g.,
$b=1$
and
$c=1$
means that osteoblasts reside “in saturation” within healthy bone. Our system of equations is
\begin{align} 0 &= \mathrm{div}(\mathbb{C}(\rho ,\sigma ,b)\varepsilon (u)) & \text{(mechanical equilibrium)}\kern30pt \end{align}
\begin{align} d_ta_i &= \mathrm{div}(D(\rho)\nabla a_i) + k_{2,i}S(\varepsilon (u)) c - k_{3,i}a_i &\text{(diffusion, generation, and}\kern26pt\nonumber\\ &&\text{decay of } i=1,2\ \text{bio-molecules)}\kern3pt \end{align}
\begin{align} d_tc &= k_6a_1a_2(1+k_7c)\bigg (1 - \frac {c}{1-\rho } \bigg) & \text{(osteoblast generation)}\kern40pt \end{align}
\begin{align} d_tb &= k_4a_1c\bigg (1-\frac {b}{1-\rho } \bigg) & \text{(bone regeneration driven by}\ a_1\text{).} \end{align}
In the above system,
$k_1,k_{2,i}, k_{3,i}, k_4, k_6, k_7\geq 0$
,
$i=1,2$
are constants that need to be determined from experiments; compare to Section 4, where we discuss certain choices. The functional relationships
$\mathbb{C}, D(\rho)$
and
$S({\cdot})$
are all required to satisfy certain technical assumptions that guarantee the well-posedness of the above system. We discuss this in detail in Section 2. For concrete examples of the functional relationships, we refer to [Reference Dondl, Poh and Zeinhofer13].
Finally, we need to specify boundary conditions. For the elastic equilibrium equation, we allow mixed boundary conditions, including the limiting cases of a pure displacement boundary condition and a pure stress boundary condition. As for the bio-active molecules, we assume that these are in saturation, i.e.,
$a_i(t,x) = 1$
,
$i=1,2$
adjacent to bone and on the rest of the boundary of
$\Omega$
we assume no-flux boundary conditions. For the initial time-point we propose
$a_i(0,x) = a_{i,0} = 0$
, for
$i=1,2$
inside of
$\Omega$
. This choice reflects the scenario of a scaffold that is not preseeded with exogenous growth factors. However, different choices of
$a_{i,0}$
are admissible and allow the model to cover, e.g., pre-seeding with osteoinductive factors. Finally, at the initial time we assume that no osteoblasts and no regenerated bone are present inside the domain of computation. In formulas, it holds for all
$i = 1,2$
\begin{align} a_i(0,x) &= 0 &\text{for all } x\in \Omega\\[-15pt]\nonumber \end{align}
\begin{align} a_i(t,x) &= 1 &\text{for all } t\in I, x \text{ adjacent to healthy bone }\\[-15pt]\nonumber \end{align}
\begin{align} D(\rho)\nabla a_i(t,x)\cdot \eta &= 0 &\text{for all } t\in I, x \text{ not adjacent to healthy bone}\\[-15pt]\nonumber \end{align}
\begin{align} (\mathbb{C}(\rho ,\sigma ,b)\varepsilon (u(t,x)))\cdot \eta &= g_N(x) & \text{on the Neumann boundary of } \Omega\\[-15pt]\nonumber \end{align}
\begin{align} u(t,x) &= g_D(x) & \text{ on the Dirichlet boundary of } \Omega\\[-15pt]\nonumber \end{align}
\begin{align} c(0,x) = b(0,x) &= 0 & \text{ for all } x\in \Omega , \end{align}
where
$\eta$
denotes the outward-pointing normal vector field on
$\partial \Omega$
. One may also consider Robin-type boundary conditions for the diffusion equations instead of equation (7). The model allows for a time-dependent choice of the mechanical loading
$g_D$
and
$g_N$
. Due to the long regeneration time horizon of approximately
$12$
months, however, it is not expedient to resolve very short time-scales of, e.g., the mechanics of physical therapy. Instead, we consider suitably time-averaged loading conditions here.
1.3. The optimization problem
In the system (1)–(4) above, the function
$\rho$
, i.e., the scaffold’s volume fraction, is a design parameter – also called control variable – that we can control in applications. For example, a given scaffold volume fraction distribution could be additively manufactured. Given a certain control variable
$\rho$
, we denote the solution of the system (1)–(4) by
\begin{equation*} y_\rho \,:\!=\, (u_\rho , a_{1,\rho }, a_{2,\rho },c_\rho ,b_\rho) \end{equation*}
to stress the dependency on
$\rho$
. We will also use the notation
$\phi (\rho) = y_\rho$
and refer to
$\phi$
as the solution operator of the system (1)–(4). Depending on the state
$y_\rho$
, we can measure the control variable’s performance by the value of an objective function
$J$
evaluated at
$\rho$
and
$y_\rho$
. We are interested in minimizing or maximizing the objective function over the set of admissible control variables. In other words, we are interested in the optimization problem of finding
\begin{equation} \underset {\rho }{\mathrm{argmin}}J(\rho ,y_\rho) \quad \text{subjected to} \quad \rho \in P, \end{equation}
where the set
$P$
encodes, for example, that
$\rho$
takes values in the unit interval (necessary for a reasonable volume fraction). The fact that, corresponding to
$\rho$
, we consider the solution
$y_\rho$
makes this a PDE-constrained optimization problem and
$\rho \in P$
introduces box constraints on the control variable.
The concrete form of
$J$
is an engineering choice. For instance, the amount of regenerated bone at a certain time-point
$T$
in the healing process should be maximized. This results in the definition
\begin{equation*} \,\hat{\!J}(\rho) = \int _\Omega b(T)\mathrm dx \end{equation*}
Care needs to be taken with respect to the functional relationships in the system (1)–(4) when choosing the amount of regenerated bone as an objective. This requires an adequate choice of the stimulus
$S$
. If
$S$
is chosen to be the Frobenius norm, the above objective function promotes very weak scaffolds as these lead to high strains and high bone growth. A more sensible choice for
$S$
in this case is to use a filter, i.e., only strains with a certain range of magnitude lead to non-vanishing values of
$S$
.
Another alternative – which we pursue here – is to minimize the temporal minimum of the elastic modulus. In the case of a soft load for the elastic equation, the elastic modulus at a time-point
$t$
is proportional to the inverse of the elastic energy
$\mathcal{E}(t)$
, i.e.,
\begin{equation*} \mathcal{E}_\rho (t) = \frac 12\int _\Omega \mathbb{C}(\rho ,\sigma (t),b(t))\varepsilon (u(t))\,:\,\varepsilon (u(t))\mathrm dx, \end{equation*}
where
$u$
and
$b$
solve the system (1)–(4) corresponding to
$\rho$
. Consequently, the maximum of
$\mathcal{E}$
over the whole regeneration process time describes the weakest state of the bone-scaffold structure during healing as a low averaged elastic modulus of the scaffold corresponds to low structural integrity. This gives rise to the objective function
\begin{equation*} \,\hat{\!J} (\rho) = \max _{t\in I}\mathcal{E}_\rho (t) \end{equation*}
and the minimization problem of finding
\begin{equation*} \rho ^* \in \underset {\rho \in P}{\mathrm{argmin}}\,\hat{\!J}(\rho). \end{equation*}
The set
$P$
encodes pointwise constraints on
$\rho$
, i.e., the necessity of enforcing
$\rho (x) \in [0,1]$
in order to be a meaningful volume fraction. The notation
$\,\hat{\!J}$
instead of
$J$
is chosen to indicate that the variables
$u$
&
$b$
appearing in the definition of
$\mathcal{E}_\rho$
are solving the system (1)–(4). Usually,
$\,\hat{\!J}$
is called the reduced objective function to distinguish it from the objective function
$J$
that does not require
$u$
and
$b$
to solve the PDE system. If we use a hard load instead of a soft load, the elastic modulus is proportional to the elastic energy; hence, the objective function becomes
\begin{equation*} \,\hat{\!J}(\rho) = \min _{t\in I}\mathcal{E}_\rho (t) \end{equation*}
and the optimization consists of finding
\begin{equation*} \rho ^* \in \underset {\rho \in P}{\mathrm{argmax}}\,\hat{\!J}(\rho), \end{equation*}
i.e., is a maximization problem.
Remark 1.
The proposed objective functions are not smooth as they involve minimizing or maximizing over
$t\in I$
. For a numerical implementation, one might therefore approximate the minimum or maximum functional by an
$L^p(I)$
norm with large value for
$-p$
or
$p$
, respectively.
2. Mathematical formulation
2.1. Notation and preliminaries
By
$X$
, we usually denote a generic Banach space and with
$X^*$
we denote its dual space. The dual pairing for
$f\in X^*$
and
$x\in X$
is denoted by
$\langle f,x \rangle _X$
. If a sequence
$(x_n)_{n\in \mathbb{N}}$
converges weakly in
$X$
to
$x$
we write
$x_n\rightharpoonup x$
. For an interval
$I=[0,T]$
and
$p\in [1,\infty ]$
we denote the Bochner space of
$p$
-integrable functions with values in
$X$
by
$L^p(I,X)$
, see for instance [Reference Diestel and Uhl11]. Further, we write
$W^{1,p}(I,X,Y)$
for the vector valued Sobolev space consisting of functions
$u\in L^p(I,X)$
with distributional derivative
$d_tu \in L^p(I,Y)$
where
$X$
and
$Y$
are Banach spaces with
$X\hookrightarrow Y$
, see for instance [Reference Boyer and Fabrie5].
We say a bounded, open set
$\Omega \subset \mathbb{R}^d$
is a Lipschitz domain if
$\overline {\Omega }$
is a Lipschitz manifold with boundary, see [Reference Grisvard17, Definition 1.2.1.2]. In the following we will denote the cube
$[{-}1,1]^n\subset \mathbb{R}^d$
by
$Q$
, its half
$\{x\in Q\mid x_d \lt 0\}$
by
$Q_-$
, the hyperplane
$\{ x\in Q \mid x_d = 0 \}$
by
$\Sigma$
and
$\{x\in \Sigma \mid x_{d-1}\lt 0\}$
by
$\Sigma _0$
. The following definition is due to Gröger, see [Reference Gröger18].
Definition 2 (Gröger regular sets). Let
$\Omega \subset \mathbb{R}^d$
be bounded and open and
$\Gamma \subset \partial \Omega$
a relatively open set. We call
$\Omega \cup \Gamma$
Gröger regular, if for every
$x\in \partial \Omega$
there are open sets
$U,V\subset \mathbb{R}^d$
with
$x\in U$
, and a bijective, bi-Lipschitz map
$\phi \,:\,U\to V$
, such that
$\phi (x) = 0$
and
$\phi (U\cap (\Omega \cup \Gamma))$
is either
$Q_-$
,
$Q_-\cup \Sigma$
or
$Q_-\cup \Sigma _0$
.
It can easily be seen that a Gröger regular set
$\Omega$
(no matter the choice
$\Gamma \subset \partial \Omega$
) is a Lipschitz domain, see [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21, Theorem 5.1]. The requirement of Gröger regularity is very mild and all applications we have in mind fall in this category. This claim is justified by the following useful characterization of Gröger regular sets in two and three dimensions that allows us to check Gröger regularity almost “by appearance”. The results are due to [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21].
Theorem 3 (Gröger regular sets in 2D, Theorem 5.2 in [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21]). Let
$\Omega \subset \mathbb{R}^2$
be a Lipschitz domain and
$\Gamma \subset \partial \Omega$
be relatively open. Then
$\Omega \cup \Gamma$
is Gröger regular if and only if
$\overline {\Gamma }\cap (\partial \Omega \setminus \Gamma)$
is finite and no connected component of
$\partial \Omega \setminus \Gamma$
consists of a single point.
Theorem 4 (Gröger regular sets in 3D, Theorem 5.4 in [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21]). Let
$\Omega \subset \mathbb{R}^3$
be a Lipschitz domain and
$\Gamma \subset \partial \Omega$
be relatively open. Then
$\Omega \cup \Gamma$
is Gröger regular if and only if the following two conditions hold
-
(i)
$\partial \Omega \setminus \Gamma$
is the closure of its interior.
-
(ii) For any
$x\in \overline {\Gamma }\cap (\partial \Omega \setminus \Gamma)$
there is an open neighbourhood
$U_x$
of
$x$
and a bi-Lipschitz map
$\phi \,:\,U_x\cap \overline {\Gamma }\cap (\partial \Omega \setminus \Gamma)\to ({-}1,1)$
.
2.2. Setting
In this section, we state the precise framework we use for the optimal control result. We begin by specifying the assumptions on the domain.
The domain. We consider a finite time interval
$I = [0,T]$
. The spatial domain
$\Omega \subset \mathbb{R}^d$
with
$d = 1,2,3$
is assumed to be an open, bounded and connected Lipschitz domain. We consider partitions of the boundary
$\partial \Omega$
, namely
\begin{gather*} \partial \Omega = \Gamma _N^e \cup \Gamma _D^e , \quad \text{and}\quad \partial \Omega = \Gamma _N^d\cup \Gamma _D^d \end{gather*}
that will be used for the elastic and the diffusion equation, respectively. Here,
$\Gamma _{N}^{e}$
and
$\Gamma _{N}^{d}$
are associated with Neumann boundary conditions, while
$\Gamma _{D}^{e}$
and
$\Gamma _{D}^{d}$
are associated with Dirichlet boundary conditions. For the partition of the elastic equation we assume
$|\Gamma _D^e| \neq 0$
. We require both
$\Omega \cup \Gamma _N^e$
and
$\Omega \cup \Gamma _N^d$
to be Gröger regular, see Definition2 or [Reference Gröger18] and [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21].
Remark 5. Note the following:
-
(i) For the elastic equation, we exclude a pure Neumann problem; however, we can include this case by passing to a suitable quotient space. We excluded this for convenience and brevity only.
-
(ii) The assumption of Gröger regularity is very mild and all desirable application settings we have in mind easily satisfy this requirement. Compare to [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21] for more information.
The control space. The set of control variables is defined to be
\begin{equation} P = \left \{ \rho \in H^2(\Omega) \mid 0 \lt c_P \leq \rho (x) \leq C_P \lt 1 \right \}\!, \end{equation}
where
$c_P$
and
$C_P$
are two fixed constants. Note that in the spatial dimensions
$d=1,2,3$
, the space
$H^2(\Omega)$
embeds into
$C^0(\Omega)$
, hence the pointwise condition imposed in the above definition is well-defined.
The state space and the equations. We denote by
$H^1_{D_e}(\Omega)$
the subspace of
$H^1(\Omega)$
with vanishing trace on
$\Gamma _D^e$
and similarly for
$H^1_{D_d}(\Omega)$
. The former is used in the context of the elasticity equation, whereas the latter is associated with the diffusion equations. Then, we consider the state space
\begin{equation*} Y = C^0\big(I,H^1_{D_e}(\Omega)\big) \times \big[H^1\big(I,H^1_{D_d}(\Omega),H^1_{D_d}(\Omega)^*\big)\cap L^4\big(I,C^0(\Omega)\big)\big]^2 \times W^{1,2}_0(I,C^0(\Omega))^2 \end{equation*}
and the space
\begin{equation*} W = L^2\big(I,H^1_{D_e}(\Omega)\big)^* \times \big[L^2\big(I,H^1_{D_d}(\Omega)\big)^*\big]^2 \times L^2(\Omega)^2 \times L^2(I,C^0(\Omega))^2. \end{equation*}
Then the state equations can be written in the form
$e(y,\rho) = 0$
with the constraint operator
\begin{equation*} e\,:\,Y\times P \to W, \quad (y,\rho) = (\tilde u,\tilde a_1,\tilde a_2, b, c, \rho)\mapsto e(y,\rho) \end{equation*}
given by
\begin{align} e(y,\rho) = \begin{pmatrix} \iint \mathbb{C}(\rho ,\sigma , b)\varepsilon (\tilde u+u_D)\,:\,\varepsilon ({\cdot})\mathrm dx\mathrm dt - \int _I \int _{\Gamma _N^e } g_N\cdot \mathrm ds\mathrm dt \\[8pt] \int _I\langle d_t\tilde a_1,\cdot \rangle _{H^1_{D_d}(\Omega)}\mathrm dt + \iint D(\rho)\nabla \tilde a_1\nabla \cdot + k_{3,1}(\tilde a_1+1)\cdot \mathrm dx\mathrm dt - \iint k_{2,1}S(\varepsilon (\tilde u +u_D)) c\cdot \mathrm dx\mathrm dt \\[8pt] \int _I\langle d_t\tilde a_2,\cdot \rangle _{H^1_{D_d}(\Omega)}\mathrm dt + \iint D(\rho)\nabla \tilde a_2\nabla \cdot + k_{3,2}(\tilde a_2+1)\cdot \mathrm dx\mathrm dt - \iint k_{2,2}S(\varepsilon (\tilde u +u_D)) c\cdot \mathrm dx\mathrm dt \\[8pt] \tilde a_1(0) + 1 \\[8pt] \tilde a_2(0) + 1 \\[5pt] d_tc - k_6(\tilde a_1 + 1)(\tilde a_2 + 1)(1+k_7c)\left (1 - \frac {c}{1 - \rho } \right) \\[8pt] d_tb - k_4(\tilde a_1 + 1)c\left (1 - \frac {b}{1 - \rho } \right) \end{pmatrix} \end{align}
Remark 6. We use the notation
\begin{equation*} a_i = \tilde a_i + 1, \quad \text{and} \quad u = \tilde u + u_D. \end{equation*}
Moreover, we frequently say a tuple
$(u, a_1, a_2, b, c, \rho)$
solves (13) for
$u\in C^0(I, H^1(\Omega))$
and
$a_i,\in H^1(I, H^1(\Omega), H^1_{D_d}(\Omega)^*)$
as an abbreviation for saying
$u=\tilde u + u_D$
and
$a_i = \tilde a_i + 1$
, where
$\tilde u$
and
$\tilde a_i$
have homogeneous boundary conditions and are the actual functions that solve (13).
Functional relationships. To make full sense of the above definition of
$e$
we still need to clarify the assumptions made on the data and functional relationships. As mentioned earlier, the function
$\sigma$
is chosen as an exponential decay
\begin{equation} \sigma (t) = \exp ({-}k_1 t). \end{equation}
We note in passing that the analysis supports general smooth
$\sigma$
that are bounded away from zero. For the material properties
$\mathbb{C}$
of the elastic equation, we require that it is a map
\begin{equation*} \mathbb{C}\,:\,\mathrm{dom}(\mathbb{C})\subset C^0(I\times \Omega)\times C^0(\Omega) \to C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s))), \quad \text{with}\quad (b,\rho)\mapsto \mathbb{C}(\rho ,\sigma ,b) , \end{equation*}
where
$\mathcal L(\mathcal M_s)$
denotes the vector space of endomorphisms of the symmetric matrices
$\mathcal M_s$
. The concrete definition of
$\mathrm{dom}(\mathbb{C})$
is not so important; however, as a minimal requirement, it should hold
\begin{equation*} \bigcup _{\rho \in P}\{ b\in C^0(I\times \Omega) \mid 0 \leq b(t,x) \leq 1 - \rho (x) \}\times \{\rho \} \subset \mathrm{dom}(\mathbb{C}). \end{equation*}
Furthermore, we need
$\mathbb{C}(\cdot ,\sigma ,\rho)$
to be Lipschitz continuous with Lipschitz constant independent of
$\rho$
and
$\mathbb{C}$
is assumed to be continuous on all of
$\mathrm{dom}(\mathbb{C})$
. Finally, we require
\begin{equation} \sup _{(b,\rho)\in \mathrm{dom}(\mathbb{C})}\lVert \mathbb{C}(\rho ,\sigma ,b) \rVert _{ C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s))) } \lt \infty \end{equation}
and
\begin{equation} \inf _{(b,\rho)\in \mathrm{dom}(\mathbb{C})}\Big[\inf _{M\in \mathcal{M}_s\setminus \{0\}} \mathbb{C}(\rho ,\sigma ,b)(t,x) M\,:\,M\Big] \geq c_{\mathbb{C}}|M|^2, \end{equation}
for a constant
$c_{\mathbb{C}}\gt 0$
uniformly in
$t\in I$
and
$x\in \Omega$
. We need a further regularity property of
$\mathbb{C}$
. We assume that
$b(t)\in C^\alpha (\Omega)$
,
$\rho \in C^\alpha (\Omega)$
for an
$\alpha \in (0,1)$
implies that the coefficient functions
\begin{equation*} C_{ijkl}(t)\,:\!=\, \left [ \mathbb{C}(\rho ,\sigma ,b) (t) \right ]_{ijkl} \end{equation*}
belong to
$C^\alpha (\Omega)$
and that there exists a constant
$C\gt 0$
not depending on
$b$
and
$\rho$
such that
\begin{equation} \lVert C_{ijkl}(t) \rVert _{C^\alpha (\Omega)} \leq C\lVert b(t) \rVert _{C^\alpha (\Omega)}\lVert \rho \rVert _{C^\alpha (\Omega)}. \end{equation}
For the boundary data
$u_D$
and
$g_N$
of the elliptic equation, we assume that
\begin{equation} g_N \in C^0(I,L^2(\partial \Omega)) \end{equation}
and that the Dirichlet boundary data are given through a function
\begin{equation} u_D \in C^0(I,H^{1+\theta }(\Omega)), \end{equation}
meaning that the boundary information can be lifted to all of
$\Omega$
such that the lift has the above regularity in time and space, where
$\theta \gt 0$
can be arbitrarily small. In practice, this is easy to verify as we mainly work with Dirichlet boundary conditions that do not vary in time. The material properties
$D(\rho)$
used in the diffusion equation are a map
\begin{equation*} D\,:\,\mathrm{dom}(D)\subset C^0(\Omega) \to L^\infty (\Omega ,\mathcal{M}_s), \quad \text{with}\quad \rho \mapsto D(\rho) \end{equation*}
that we require to be continuous with respect to the uniform norm on
$\mathrm{dom}(D)$
. The domain of
$D$
will usually satisfy
\begin{equation*} \{ \rho \in C^0(\Omega) \mid 0 \lt c_P \leq \rho (x) \leq C_P \lt 1 \} \subset \mathrm{dom}(D), \end{equation*}
where
$c_P$
and
$C_P$
are the positive constants appearing in the definition of
$P$
. We also require
$D$
to be uniformly elliptic independently of
$\rho \in \mathrm{dom}(D)$
, i.e.,
\begin{equation} \inf _{\rho \in \mathrm{dom}(D)}\Big[ \inf _{\xi \in \mathbb{R}^d\setminus \{0\}} D(\rho)\xi \cdot \xi \Big] \geq c_D |\xi |^2, \end{equation}
for a constant
$c_D \gt 0$
. Finally, for the function
$S({\cdot})$
, we assume that it is given through a map on matrices
\begin{equation*} S({\cdot})\,:\,\mathbb{R}^{d\times d} \to [0,\infty) \end{equation*}
that we require to be Lipschitz and to obey an estimate of the form
\begin{equation} S(A) \leq C_1|A| + C_2\quad \text{for all }A\in \mathbb{R}^{d\times d}, \end{equation}
where
$C_1, C_2 \gt 0$
and
$|A|$
denotes the Euclidean (or any) norm of a matrix. Furthermore, we need
$S$
to be continuous, more precisely, we assume that if
$(v_k)\subset L^2(\Omega ,\mathbb{R}^{d\times d})$
is a sequence, then it holds
\begin{equation} v_k \to v \quad \text{in }L^2(\Omega ,\mathbb{R}^{d\times d}) \quad \Rightarrow \quad S(v_k) \to S(v) \quad \text{in }L^2(\Omega). \end{equation}
We state the main result concerning the well-posedness of the PDE-ODE system.
Theorem 7.
Assume that the setting described in this section holds. Then, for every
$\rho \in P$
there exists a unique solution
$y=(\tilde u, \tilde a_1, \tilde a_2, b, c)\in Y$
satisfying
$e(y,\rho) = 0$
, i.e., solving the state equations (13).
Proof. This follows almost as an application of Theorem 3.2 in [Reference Dondl, Poh and Zeinhofer13]. Note that our assumptions here are slightly stronger, so the requirements in [Reference Dondl, Poh and Zeinhofer13] are trivially satisfied. The only extension to the results in [Reference Dondl, Poh and Zeinhofer13] is to show the improved integrability for the functions
$a_1$
and
$a_2 \in L^4(I,C^0(\Omega))$
. To this end, we use the estimate
$(3.11)$
in [Reference Dondl, Poh and Zeinhofer13], which shows that the right-hand sides of the diffusion equations satisfy
\begin{equation*} k_{2,i}S(\varepsilon (\tilde u +u_D)) c \in L^\infty (I,L^2(\Omega)). \end{equation*}
This allows to apply the maximal
$L^p$
regularity result in Lemma 21 and obtain
$\tilde a_i \in L^p(I,C^0(\Omega))$
for all
$p\in [2,\infty)$
.
2.3. Objective function
Here we formulate the class of objective functions we are able to treat in the setting of the optimal control result. For every time-point
$t\in I$
and state control pair
$(y,\rho) \in Y \times P$
, we consider the elastic energy
\begin{equation} \mathcal{E}\,:\,Y\times P \to C^0(I), \quad \mathcal{E}(y,\rho)(t) = t\mapsto \frac 12 \int _\Omega \mathbb{C}(\rho ,\sigma (t),b(t)) \varepsilon (u(t))\,:\,\varepsilon (u(t))\mathrm dx. \end{equation}
For most of our objective functions, we desire
$\mathcal{E}$
to take values in
$C^0(I)$
, as we want to have access to point evaluations. This is the reason to require the continuity of the solutions to the elastic equation in the definition of
$Y$
. Primarily, we are interested in the reduced elastic energy
$\hat {\mathcal{E}}$
, that is, we are interested in
$\mathcal{E}(y,\rho)$
only when
$(y,\rho)$
solves the system of equations, i.e., when it holds
$e(y,\rho) = 0$
. We define
\begin{equation} \hat {\mathcal{E}}\,:\, P \to C^0(I), \quad \hat {\mathcal{E}}(\rho) = \mathcal{E}(y_\rho ,\rho) \end{equation}
and here it holds
$e(y_\rho ,\rho) = 0$
. We provide now the proof that
$\mathcal{E}$
takes values in
$C^0(I)$
.
Lemma 8.
For all
$(y,\rho)\in Y\times P$
, we have
$\mathcal{E}(y,\rho)\in C^0(I)$
. If it holds
$e(y,\rho) = 0$
and
$\tilde u(t) +u_D(t) \neq 0$
, then
$\mathcal{E}(y,\rho)(t) \gt 0$
.
Proof. As
$\tilde u+u_D \in C^0(I, H^1(\Omega))$
, and the material tensor
$\mathbb{C}$
is in the space
$C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))$
, it follows that
\begin{equation*} \mathbb{C}(\rho ,\sigma ,b) \varepsilon (\tilde u+u_D)\,:\,\varepsilon (\tilde u+u_D) \in C^0(I,L^1(\Omega)). \end{equation*}
Using the
$L^1(I)$
continuity of integration, we get
$\mathcal{E}(y,\rho)\in C^0(I)$
. Now, let
$e(y,\rho) = 0$
. We can estimate
\begin{equation*} \mathcal{E}(y,\rho)(t) \geq c \lVert \tilde u(t) + u_D \rVert ^2_{H^1(\Omega)}, \end{equation*}
with the constant
$c \gt 0$
depending on the constant appearing in Korn’s inequality and the ellipticity constant
$c_{\mathbb{C}}$
. As it holds
$e(y,\rho) = 0$
, for every
$t\in I$
the function
$\tilde u(t) + u_D(t)$
solves an elastic equation, hence can only vanish if the boundary conditions are homogeneous for this time-point, which leads to
$\tilde u(t) + u_D(t) = 0$
. This is excluded in the statement of the lemma and the proof is complete.
We state now the structural assumption we impose for our admissible objective functions.
Assumption 9.
Let
$\mathcal{F}\,:\,\mathrm{dom}(\mathcal{F})\subset C^0(I) \to \mathbb{R}$
be a continuous map and assume that the domain of
$\mathcal{F}$
satisfies
\begin{equation} \{ v \in C^0(I) \mid v(t) \gt 0 \text{ for all }t\in I \} \subset \mathrm{dom}(\mathcal{F}). \end{equation}
Furthermore, let
$\mathcal{G}\,:\,C^0(I\times \Omega) \to \mathbb{R}$
be a continuous function. Using the elastic energy
$\mathcal{E}$
and functionals
$\mathcal{F}$
,
$\mathcal{G}$
as above, we define the prototypical objective function as
\begin{equation*} J\,:\, Y \times P \to \mathbb{R}, \quad J(y,\rho) = \mathcal{F}\left (\mathcal{E}(y,\rho) \right) + \mathcal{G}(b) \end{equation*}
in case the domain of
$\mathcal{F}$
allows
$\mathcal{E}(y,\rho)$
as an argument. The function
$b$
denotes the bone component of the state variable
$y$
. More importantly, we define the reduced objective
\begin{equation*} \,\hat{\!J} \,:\, P \to \mathbb{R}, \quad \,\hat{\!J}(\rho) = \mathcal{F}\left (\mathcal{E}(\phi (\rho),\rho) \right) + \mathcal{G}(b), \end{equation*}
where
$\phi \,:\,P \to Y$
denotes the control-to-state operator, i.e., for a given
$\rho$
,
$\phi (\rho)$
solves the associated system. Note that the assumption (25) together with Lemma 8 guarantees that
$\mathcal{E}(\phi (\rho),\rho)$
is an admissible argument of
$\mathcal{F}$
. Finally, we assume that
$\,\hat{\!J}$
is bounded from below if we are interested in a minimization problem and we assume
$\,\hat{\!J}$
to be bounded from above if we are interested in maximization.
Remark 10. We discuss how the examples discussed in Section 1.3 fall in the abstract setting described above.
-
(i) Choosing the minimum (or maximum) functional
for
\begin{equation*} \min \,:\, C^0(I) \to \mathbb{R}, \quad v\mapsto \min _{t\in I}v(t) \end{equation*}
$\mathcal{F}$
is conforming with Assumption 9 as clearly
$\mathrm{min}$
and
$\mathrm{max}$
are continuous functionals on
$C^0(I)$
.
-
(ii) Smooth approximations of the minimum and the maximum are given by
$L^p(I)$
norms with large values of
$|p|$
. A positive value for
$p$
serves as an approximation of the maximum and a negative value is suitable for the approximation of the minimum. In the latter case, i.e.,
$p\lt 0$
, one chooses
It is straightforward to show that
\begin{equation*} \mathrm{dom}\left (\lVert \cdot \rVert _{L^p(I)} \right) :\!=\{ v \in C^0(I) \mid v(t) \gt 0 \text{ for all }t\in I \}. \end{equation*}
$\lVert \cdot \rVert _{L^p(I)}$
is continuous with respect to the uniform norm, also for negative exponents. In fact, it is even Fréchet differentiable.
-
(iii) The choice
corresponds to the objective of regenerated bone at time
\begin{equation*} \mathcal{G}(b) = \int _\Omega b(T) \mathrm dx, \end{equation*}
$T$
. Clearly,
$\mathcal{G}$
is continuous and evaluating
$\mathcal{G}$
only at functions
$b$
that solve the state equations shows that
$\mathcal{G}$
is bounded.
3. Main results
Our main result establishes the existence of an optimal control in the set
$P\subset H^2(\Omega)$
given the objective function
$\,\hat{\!J}$
is regularized by an
$H^2(\Omega)$
norm.
Theorem 11 (Optimal control). Assume we are in Setting 2.2 and let
$\eta \gt 0$
be fixed. Then there exists a minimizer
$\rho ^* = \rho ^*(\eta) \in P$
to the regularized objective
\begin{equation*} \,\hat{\!J}(\rho ^*) + \eta \lVert \rho ^* \rVert _{H^2(\Omega)} = \inf _{\rho \in P}\big[ \,\hat{\!J}(\rho) + \eta \lVert \rho \rVert ^2_{H^2(\Omega)} \big]. \end{equation*}
In order to incorporate the pointwise constraint encoded in the definition of the control space
$P$
, see (12), in a numerical simulation one can use a soft penalization. This usually corresponds to a continuous functional
$\mathcal{K}\,:\,C^0(\Omega)\to [0,\infty)$
. Also, in this setting, we can establish the existence of an optimal control.
Corollary 12.
Assume we are in Setting 2.2 and let
$\mathcal{K}\,:\,C^0(\Omega)\to [0,\infty)$
be a continuous, non-negative functional. Then there exists an optimal control
$\rho ^\dagger = \rho ^\dagger (\eta ,\mathcal{K}) \in H^2(\Omega)$
to the regularized and penalized objective, i.e.,
\begin{equation*} \,\hat{\!J}(\rho ^\dagger) + \eta \lVert \rho ^\dagger \rVert ^2_{H^2(\Omega)} + \mathcal{K}(\rho ^\dagger) = \inf _{\rho \in H^2(\Omega)}\big[ \,\hat{\!J}(\rho) + \eta \lVert \rho \rVert ^2_{H^2(\Omega)} + \mathcal{K}(\rho) \big]. \end{equation*}
Remark 13. A few comments regarding the above results are in order.
-
(i) For some objectives we might be interested in a maximizer rather than a minimizer. In this case, one subtracts the regularizer
$\eta \lVert \cdot \rVert _{H^2(\Omega)}$
and the soft penalty
$\mathcal{K}$
and the results are still valid. For brevity, we discuss only minimization problems in the remainder.
-
(ii) As discussed in Section 2.3 , we have some freedom in the choice of
$\,\hat{\!J}$
. From a modelling perspective, a maximum or minimum over all time-points of the elastic energy seems reasonable. On the other hand, for the numerical treatment a smooth approximation thereof is preferable, e.g., an
$L^p(I)$
norm. Note that all these choices are covered by our main result.
-
(iii) The Tikhonov penalization term
$\eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
is artificial. It serves to generate compactness of minimizing sequences and an optimal control result without this term seems out of reach.
-
(iv) It is presently unclear to us if the optimal control problem possesses a unique solution.
The strategy to prove Theorem 11 and Corollary 12 is the direct method of the calculus of variations and crucially relies on rather specific regularity properties of the diffusion equations and the elastic equation that imply convenient compact embeddings. The technical results concerning these regularity properties are established in Appendix A. In this section, we assume the implications of the compact embeddings and show how this leads to a proof of Theorem 11. We stress that the mixed boundary conditions, rough coefficients and jump initial conditions are responsible for the technical difficulties.
Proposition 14.
Assume we are in Setting 2.2. Let
$(\rho _k)\subset P$
be a minimizing sequence for
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
and denote by
$(u_k) \subset C^0(I,H^1(\Omega))$
,
$(a_{1, k}),(a_{2, k}) \subset H^1(I,H^1(\Omega),H^1_{D_d}(\Omega)^*) \cap L^4(I, C^0(\Omega))$
and
$(b_k), (c_k) \subset W^{1,2}(I,C^0(\Omega))$
the corresponding solutions to the system 13, refer also to Remark 6
. Assume that there is a common subsequence (not relabelled) of
$(\rho _k), (u_k), (a_{1, k}), (a_{2, k}), (b_k), (c_k)$
and elements
$\rho ^* \in P$
,
$u^* \in C^0(I,H^1(\Omega))$
,
$a_1^*$
,
$a_2^* \in H^1(I,H^1(\Omega),H^1_D(\Omega)^*)$
and
$b^*, c^* \in W^{1,2}(I, C^0(\Omega))$
such that
-
(A1)
$\rho _k \to \rho ^*$
in
$C^0(\Omega)$
and
$\rho _k \rightharpoonup \rho ^*$
in
$H^2(\Omega)$
, -
(A2)
$u_k \to u^*$
in
$C^0(I, H^1(\Omega))$
, -
(A3)
$a^i_k \rightharpoonup a_i^*$
in
$H^1(I,H^1(\Omega), H^1_D(\Omega)^*)$
,
$i=1,2$
-
(A4)
$b_k \to b^*$
in
$C^0(I\times \Omega)$
-
(A5)
$c_k \to c^*$
in
$C^0(I\times \Omega)$
then
$(\rho ^*, u^*, a_1^*, a_2^*, b^*)$
solves the system 13 and
$\rho ^*$
is minimizer of
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
over the set
$P$
, i.e., satisfies
\begin{equation*} \,\hat{\!J}(\rho ^*) + \eta \lVert \rho ^* \rVert _{H^2(\Omega)} = \inf _{\rho \in P}\big[ \,\hat{\!J}(\rho) + \eta \lVert \rho \rVert ^2_{H^2(\Omega)} \big]. \end{equation*}
Proof. There are two things to show. First, we need to guarantee that the tuple
$(\rho ^*,u^*,a_1^*, a_2^*,b^*)$
still solves the system of equations (13). And secondly, we need to prove that
$\rho ^*$
is in fact a minimizer. We start with the second point, assuming for the moment that
$(\rho ^*,u^*,a_1^*, a_2^*,b^*)$
solves the correct equations. We show that it holds
\begin{equation*} \,\hat{\!J}(\rho ^*) + \eta \lVert \rho ^* \rVert _{H^2(\Omega)}^2 \leq \liminf _{k\to \infty }\big[ \,\hat{\!J}(\rho _k) + \eta \lVert \rho _k \rVert ^2_{H^2(\Omega)} \big] = \min _{\rho \in P}\big[ \,\hat{\!J}(\rho) + \eta \lVert \rho \rVert ^2_{H^2(\Omega)} \big], \end{equation*}
that is, the classical lower semi-continuity property required in the application of the direct method of the calculus of variations. Clearly, the map
\begin{equation*} H^2(\Omega)\to \mathbb{R}, \quad \rho \mapsto \eta \lVert \rho \rVert _{H^2(\Omega)}^2 \end{equation*}
is convex and norm continuous, hence weakly lower semi-continuous, that is, it holds
\begin{equation*} \eta \lVert \rho ^* \rVert _{H^2(\Omega)}^2 \leq \liminf _{k\to \infty }\eta \lVert \rho _k \rVert _{H^2(\Omega)}^2 \end{equation*}
by the assumption
$\rho _k\rightharpoonup \rho ^*$
in
$H^2(\Omega)$
on the minimizing sequence. To proceed, remember our structural assumption on the objective function, i.e.,
\begin{equation*} \,\hat{\!J} = \mathcal{F}(\hat {\mathcal{E}}(\rho)) + \mathcal{G}(b), \end{equation*}
where
$\mathcal{F}\,:\, \mathrm{dom}(\mathcal F) \subset C^0(I)\to \mathbb{R}$
and
$\mathcal{G}\,:\,C^0(I\times \Omega)\to \mathbb{R}$
are assumed to be continuous. Thus it suffices to show that
$\mathcal{E}(\rho _k)\to \mathcal{E}(\rho ^*)$
in
$C^0(I)$
. For convenience, let us now set
$\mathbb{C}^* = \mathbb{C}(\rho ^*,\sigma ,b^*)$
and
$\mathbb{C}_k = \mathbb{C}(\rho _k,\sigma ,b_k)$
. We then compute
\begin{align*} \lVert \hat {\mathcal{E}}(\rho _k) - \hat {\mathcal{E}}(\rho ^*) \rVert _{C^0(I)} &= \frac 12 \left \lVert \int _\Omega \left [\mathbb{C}_k-\mathbb{C}^*\right ]\varepsilon (u_k)\,:\,\varepsilon (u_k) + \mathbb{C}^*\varepsilon (u_k-u^*)\,:\,\varepsilon (u_k) + \mathbb{C}^*\varepsilon (u^*)\varepsilon (u_k-u^*)\mathrm dx \right \rVert _{C^0(I)} \\[2pt] &\leq \frac 12\lVert \mathbb{C}_k - \mathbb{C}^* \rVert _{C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))} \lVert \varepsilon (u_k) \rVert ^2_{C^0(I,L^2(\Omega))} \\[2pt] & \quad + \frac 12\lVert \mathbb{C}^* \rVert _{C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))} \lVert \varepsilon (u_k - u^*) \rVert _{C^0(I,L^2(\Omega))} \lVert \varepsilon (u_k) \rVert _{C^0(I,L^2(\Omega))} \\[2pt] &\quad + \frac 12\lVert \mathbb{C}^* \rVert _{C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))} \lVert \varepsilon (u^*) \rVert _{C^0(I,L^2(\Omega))} \lVert \varepsilon (u_k - u^*) \rVert _{C^0(I,L^2(\Omega))}. \end{align*}
Using the continuity assumption for
$\mathbb{C}$
and the convergence
$b_k\to b^*$
in
$C^0(I\times \Omega)$
and
$\rho _k\to \rho ^*$
in
$C^0(\Omega)$
we get that
\begin{equation*} \lVert \mathbb{C}_k - \mathbb{C}^* \rVert _{C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))} \to 0. \end{equation*}
Furthermore, the convergence
$u_k\to u^*$
in
$C^0(I,H^1(\Omega))$
implies both a bound on
$\lVert \varepsilon (u_k) \rVert$
and the convergence
\begin{equation*} \lVert \varepsilon (u_k - u^*) \rVert _{C^0(I,L^2(\Omega))}. \end{equation*}
Hence, we established
$\hat {\mathcal E}(\rho _k)\to \hat {\mathcal{E}}(\rho ^*)$
and conclude
\begin{equation*} \,\hat{\!J}(\rho ^*) + \eta \lVert \rho ^* \rVert ^2_{H^2(\Omega)} \leq \lim _{k\to \infty }\,\hat{\!J}(\rho _k) + \liminf _{k\to \infty }\eta \lVert \rho _k \rVert ^2_{H^2(\Omega)} \leq \liminf _{k\to \infty }\big[ \,\hat{\!J}(\rho _k) + \eta \lVert \rho _k \rVert _{H^2(\Omega)}^2 \big] \end{equation*}
which settles the claim.
We still need to show that
$(\rho ^*,u^*,a_1^*, a_2^*,b^*)$
is in fact a solution to the system 13. For the elastic equation, we consider an arbitrary test function
$\varphi \in L^2(I,H^1_{D_e}(\Omega))$
\begin{equation*} \iint \mathbb{C}(\rho _k,\sigma ,b_k)\varepsilon (u_k)\,:\,\varepsilon (\varphi)\mathrm dx\mathrm dt = \int _I \int _{\partial \Omega }g_N\varphi \mathrm ds\mathrm dt \end{equation*}
and the continuity assumption on
$\mathbb{C}$
and the convergence assumed for
$\rho _k$
,
$b_k$
and
$u_k$
are by far sufficient to pass to the limit.
In the same spirit, we consider the diffusion equations with a test function
$\varphi \in L^2(I,H_{D_d}(\Omega))$
\begin{equation*} \int _I\langle d_ta_{i,k},\varphi \rangle _{H^1_{D_d}(\Omega)}\mathrm dt + \iint D(\rho _k)\nabla a_{i,k}\nabla \varphi + k_3 a_{i,k} \varphi \mathrm dx\mathrm dt = \iint k_2 S(\varepsilon (u_k)) c_k\varphi \mathrm dx\mathrm dt, \quad i=1,2. \end{equation*}
For the left-hand side of the diffusion equations, we can easily pass to the limit by the weak convergence of
$a_{i,k}$
and the strong convergence of
$D(\rho _k)$
that we have available through the continuity assumption on
$D$
and
$\rho _k\to \rho ^*$
in
$C^0(\Omega)$
. For the right-hand side, we use the implication
\begin{equation*} u_k\to u^*\text{ in }C^0(I,H^1(\Omega))\quad \Rightarrow \quad S(\varepsilon (u_k)) \to S(\varepsilon (u^*)) \text{ in }L^2(\Omega). \end{equation*}
Hence, the limit for the diffusion equations can also be correctly identified. To establish the initial condition of the limit, consider the continuous linear map
\begin{equation*} H^1\big(I,H^1_{D_d}(\Omega), H^1_{D_d}(\Omega)^*\big) \to C^0(I,L^2(\Omega)) \to L^2(\Omega), \quad a\mapsto a(0). \end{equation*}
Using the weak sequential continuity of continuous linear maps shows that
$a^*(0)$
vanishes, as desired.
To pass to the limit in the cell ODE, we look at its fixed-point equation
\begin{equation*} c_k(t) = \int _0^t k_6 a_{1,k}(s)a_{2,k}(s) (1 + k_7c_k(s))\left (1 - \frac {c_k(s)}{1 - \rho _k}\right)\mathrm ds, \end{equation*}
which holds in the space
$C^0(\Omega)$
, for all
$t\in I$
. Multiplying the above equation by a smooth test function
$\varphi \in C_c^\infty (\Omega)$
and integrating over
$\Omega$
yields for the left-hand side of the above equation
\begin{equation*} \int _\Omega c_k(t)\varphi \mathrm dx \to \int _\Omega c^*(t)\varphi \mathrm dx \quad \text{with}\quad k\to \infty . \end{equation*}
The convergence
$c_k \to c^*$
in the space
$C^0(I\times \Omega)$
suffices by far for the above limit passage. Before we treat the limit of the right-hand side, we note that the compactness result of Aubin-Lions, see, for instance [Reference Simon34], provides the compact embedding
\begin{equation*} H^1\big(I,H^1_D(\Omega), H^1_D(\Omega)^*\big) \hookrightarrow \hookrightarrow L^2(I,L^2(\Omega)) \end{equation*}
which is essentially due to the fact that the space triple
$(H^1_D(\Omega),L^2(\Omega), H^1_D(\Omega)^*)$
satisfies the requirements of the Ehrling Lemma, being in turn guaranteed by the Rellich–Kondrachov compactness result that provides the compact embedding of
$H^1_D(\Omega)$
into
$L^2(\Omega)$
. Note that the boundary regularity for
$\Omega$
is chosen to support the Rellich–Kondrachov theorem. Hence we get the convergence
\begin{equation*} a_{1,k} a_{2,k} \to a^*_1 a^*_2 \quad \text{in}\quad L^1(I,L^1(\Omega))\cong L^1(I\times \Omega). \end{equation*}
Using the above convergence and the convergence of
$c_k\to c^*$
in
$C^0(I\times \Omega)$
and
$\rho _k \to \rho ^*$
in
$C^0(\Omega)$
, we compute, employing Fubini’s theorem and pass to the limit
\begin{align*} \int _\Omega \int _0^t k_6 a_{1,k} a_{2,k} (1 + k_7c_k)\left (1 - \frac {c_k}{1 - \rho _k}\right)\mathrm ds \varphi \mathrm dx &= \int _0^t\int _\Omega k_6 a_{1,k} a_{2,k} (1 + k_7c_k)\left (1 - \frac {c_k}{1 - \rho _k}\right)\varphi \mathrm dx \mathrm ds \\[3pt] &\to \int _0^t\int _\Omega k_6a^*_1 a^*_2 (1 + k_7c^*)\left (1 - \frac {c^*}{1 - \rho ^*}\right)\varphi \mathrm dx \mathrm ds \\[3pt] &= \int _\Omega \int _0^t k_6 a^*_1 a^*_2 (1 + k_7c^*)\left (1 - \frac {c^*}{1 - \rho ^*}\right) \mathrm ds \varphi \mathrm dx \end{align*}
Inferring the fundamental lemma of the calculus of variations, we obtain
\begin{equation*} c^*(t) = \int _0^t k_6 a_1^* a_2^* (1+k_7c^*)\left (1 - \frac {c^*}{1-\rho ^*} \right)\mathrm ds \end{equation*}
for every
$t \in I$
. This implies that
$c^*$
satisfies the correct limit equation. Obviously, we can repeat the same argument to guarantee that
$b^*$
satisfies an appropriate limit equation.
Remark 15.
Via discussing the requirements
$(A1){-}(A5)$
above, we give a rough idea of their proof.
-
(i) The fact that
$J$
is bounded from below implies that the regularization term
$\eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
automatically leads to an
$H^2(\Omega)$
bound on any minimizing sequence
$(\rho _k)\subset P$
. Thus there exists
$\rho ^* \in P$
and a (not re-labelled) subsequence
$(\rho _k)$
with
$\rho _k \rightharpoonup \rho ^*$
in
$H^2(\Omega)$
. Employing the compactness
that holds for three spatial dimensions, this implies the desired convergence
\begin{equation*} H^2(\Omega) \hookrightarrow \hookrightarrow C^0(\Omega) \end{equation*}
$\rho _k \to \rho ^*$
in
$C^0(\Omega)$
.
-
(ii) A uniform bound in
$C^0(I,H^1(\Omega))$
norm of the sequence
$(u_k)$
is easily established as Lemma 16 shows. However, this does not provide assumption (A2), which can only be achieved through a compactness argument. In fact – given Hölder continuous coefficients functions of
$\mathbb{C}(\rho _k,\sigma , b_k)$
– one is able to show that for every
$t\in I$
the solution
$u_k(t)$
is in
$H^{1+\theta }(\Omega)$
for a sufficiently small
$\theta \gt 0$
as an application of the main theorem of [Reference Haller-Dintelmann, Meinlschmidt and Wollner20]. Compare also to Lemma 18 for a discussion of the applicability of this result. Then, given the relative compactness of the sequences
$(b_k)$
in
$C^0(I\times \Omega)$
and
$(\rho _k) \subset C^0(\Omega)$
, one can apply a vector-valued version of the Arzelà-Ascoli theorem to derive the relative compactness of
$(u_k)$
in
$C^0(I,H^1(\Omega))$
. As discussed in (iv), the compactness of
$(b_k)$
relies on a Hölder regularity result for diffusion equations.
-
(iii) Similarly, a uniform bound for the sequences
$(a_{i,k})$
in
$H^1(I,H^1(\Omega),H^1_D(\Omega)^*)$
norm can be established by standard computations, thus implying the desired existence of
$a_i^*$
and corresponding subsequence. We provide the details in Lemma 20
. -
(iv) The existence of a subsequence
$(b_k)$
and
$b^*\in C^0(I\times \Omega)$
with
$b_k\to b^*$
in
$C^0(I\times \Omega)$
requires the biggest effort. We achieve this by deriving a
$W^{1,2}(I,C^\alpha (\Omega))$
bound on
$(b_k)$
for an
$\alpha \in (0,1)$
. Investigating the structure of the cell and bone ODEs, we see that such a regularity and bound can only be established if we are able to show that the sequences
$(a_{i,k})$
are bounded in
$L^2(I,C^\alpha (\Omega))$
. It is this regularity and boundedness result for the diffusion equation on which the whole proof rests, we state it in Lemma 21
, but the derivation of this result is the topic of [Reference Dondl and Zeinhofer14].Coming back to the boundedness of
$(b_k)$
in
$W^{1,2}(I,C^\alpha (\Omega))$
, note that this implies the desired existence of
$b^*\in C^0(I\times \Omega)$
together with a subsequence
$b_k \to b^*$
in
$C^0(I\times \Omega)$
via the embeddings
\begin{equation*} W^{1,2}(I,C^\alpha (\Omega)) \hookrightarrow C^\beta (I,C^\alpha (\Omega)) \hookrightarrow C^{\min (\alpha ,\beta)}(I\times \Omega)\hookrightarrow \hookrightarrow C^0(I\times \Omega). \end{equation*}
-
(vi) To summarize:
$(A1)$
is clear,
$(A3)$
is established in Lemma 20
,
$(A2)$
,
$(A4)$
and
$(A5)$
rely on the regularity result for diffusion equations stated in Lemma 21 and the main result of [Reference Haller-Dintelmann, Meinlschmidt and Wollner20]. The derivation of the
$W^{1,2}(I,C^\alpha (\Omega))$
bound for
$(b_k)$
is carried out in Lemma 25
, the bound for
$(c_k)$
in Lemma 24
.
4. Simulations
In this section, we present numerical simulations of optimal scaffold density distributions. Our main application example is large tibial defects and we are especially interested in stress shielding effects caused by external fixation of the scaffold. Our numerical findings indicate that a three-dimensional scaffold density optimization is of substantial importance in the mitigation of such stress shielding effects, which usually cause reduced bone regeneration in the vicinity of the external fixateur, see, e.g., [Reference Viateau, Guillemin, Bousson, Oudina, Hannouche, Sedel, Logeart‐Avramoglou and Petite38, Figure 5c].
4.1. Stress shielding
Bone adapts according to the mechanical environment it is subjected to. This important property of bone is well known and commonly referred to as Wolff’s law, see [Reference Wolff40]. It has far ranging consequences for bone tissue engineering. More precisely, prosthetic implants are often made of stiffer materials than bone and thus change the mechanical environment in their vicinity. This often leads to bone regions that are subjected to less stress and consequently bone resorption when compared to a healthy bone, a phenomenon known as stress shielding, which has been extensively studied, e.g., in the context of total hip arthroplasty, see [Reference Arabnejad, Johnston, Tanzer and Pasini2, Reference Behrens, Wirth, Windhagen, Nolte, Meyer-Lindenberg and Bouguecha3, Reference Huiskes, Weinans and Van Rietbergen23, Reference Sumner and Galante36]. The bone resorption in the vicinity of the prosthetic implant can lead to serious complications such as periprosthetic fracture and aseptic loosening and revision surgeries – if so needed – can be complicated, we refer to [Reference Arabnejad, Johnston, Tanzer and Pasini2].
It is to be expected that stress shielding effects also play an important role in scaffold-mediated bone growth, for example caused through the external fixation of the scaffold by a metal plate. This leads to under-loading in the vicinity of the fixating element. To be able to quantify these effects, it is crucial to use a three-dimensional computational model; a one-dimensional simplification, as for instance discussed by [Reference Poh, Valainis, Bhattacharya, van Griensven and Dondl32], cannot resolve the asymmetries that induce the effect.
4.2. The computational model
Our concrete model setup is identical to the one presented in [Reference Dondl, Poh and Zeinhofer13] as far as the state equations are concerned. For the reader’s convenience, we briefly repeat the state equations and boundary conditions
\begin{align} 0 &= \mathrm{div}(\mathbb{C}(\rho ,\sigma , b)\varepsilon (u)) &{\mathbb{C}(\rho ,\sigma , b) = b\mathcal{C}_b + \rho \sigma \mathcal{C}_\rho , \ \sigma (t)=\exp ({-}0.1t)}\\[-13pt]\nonumber \end{align}
\begin{align} d_ta_1 &= \mathrm{div}(k_5(1-\rho) \nabla a_1) + k_{2,1}|\varepsilon (u)|c - k_{3,1}a_1 &{k_5=260, \ k_{2,1}=10500, \ k_{3,1}=16.}\kern46pt\\[-13pt]\nonumber \end{align}
\begin{align} d_ta_2 &= \mathrm{div}(k_5(1-\rho) \nabla a_2) + k_{2,2}|\varepsilon (u)|c - k_{3,2}a_2 &{k_{2,2}=5250, \ k_{3,1}=8.}\kern98pt \end{align}
\begin{align} d_tc &= k_6a_1a_2(1 + k_7c)\bigg (1 - \frac {c}{1 - \rho } \bigg)\kern42pt &{k_6=0.39, \ k_7=0.07.}\kern98pt\\[-13pt]\nonumber \end{align}
\begin{align} d_tb &= k_4a_1c\bigg (1 - \frac {b}{1 - \rho } \bigg) &{k_4=0.14}\kern145pt \end{align}
where
$\mathcal{C}_b$
and
$\mathcal{C}_\rho$
are the material tensors for PCL and bone modelled as linear isotropic materials. For the Lamé constants, we use
\begin{gather*} (\lambda _\rho ,\mu _\rho)=(1.97,0.17)\mathrm{GPa}, \quad (\lambda _b,\mu _b)=(2.88,1.92)\mathrm{GPa}, \end{gather*}
for PCL and bone, respectively. We use the same boundary conditions as in [Reference Dondl, Poh and Zeinhofer13] with the exception of the elastic equation, which is now subject to pure stress boundary conditions with a constant surface traction stemming from a force of
$0.3\mathrm{kN}$
which is applied to the top and bottom of the cylindrical domain. We propose to view this as a typical force occurring during the regeneration time. The bioactive molecules
$a_1,a_2$
are assumed to be in saturation adjacent to the initial, healthy bone matrix at the top and bottom of the domain and a scenario without preseeding throughout the domain (i.e., a zero initial condition) is considered. For a detailed reasoning for the choices of parameters, boundary conditions and functional relationships, we refer to [Reference Dondl, Poh and Zeinhofer13].
As an objective function to measure scaffold performance, we use the maximum over the temporal evolution of the scaffold-bone composite’s elastic energy. Due to the stress boundary condition in the numerical experiments, the reciprocal of the elastic energy is proportional to the elastic modulus of the scaffold-bone system: a reasonable measure of its mechanical integrity. The optimization goal is to minimize this temporal maximum while respecting the state equations and an additional constraint on
$\rho$
to not take values outside the unit intervalFootnote
1
. In formulas, we denote by
$\mathcal{E}$
the elastic energy
\begin{equation*} \mathcal{E}(y,\rho)(t) = \frac 12\int _\Omega \mathbb{C}(\rho (x),\sigma (t),b(t,x))\varepsilon (u(t,x))\,:\,\varepsilon (u(t,x))\mathrm{d}x \end{equation*}
where
$y=(u,a_1,a_2,c,b)$
is the state variable. The minimization problem is the task to find
\begin{equation} \rho \in \mathrm{argmin}\big[ \max _{t\in I}\mathcal{E}(y,\rho)(t) \big], \quad \text{subjected to }e(y,\rho)=0\quad \text{and}\quad \rho \in P, \end{equation}
where
$P$
encodes that
$\rho$
is bounded away from zero and one. Numerically, we replace the temporal maximum by an
$L^p(I)$
norm with
$p=5$
to smoothly approximate it. The pointwise constraint
$\rho \in P$
is treated by a soft penalty and
$e(y,\rho)=0$
by the adjoint method.
4.3. Numerical implementation
The numerical realization of the PDE constrained optimization problem is based on the adjoint approach, see for instance [Reference Hinze, Pinnau, Ulbrich and Ulbrich22] for a derivation of the method. This means that the constraint
$e(y,\rho)=0$
(in the notation of Section 2) is parametrized by the solution operator
$\rho \mapsto \phi (\rho)$
satisfying
$e(\phi (\rho),\rho)=0$
eliminating the constraint in the optimization. Computing the derivative of the reduced objective with respect to
$\rho$
yields an adjoint equation that is structurally similar to the state equations (13). Having access to the derivative of the reduced objective, we use an
$L^2(\Omega)$
gradient flow in order to solve the optimization problem (31). In the numerical implementation, we do not employ the
$H^2$
regularization of the objective function
$\,\hat{\!J}$
. As this leads to reasonable results, more sophisticated optimization algorithms were not deemed necessary.
We use the Computational Geometry Algorithms Library ([Reference Boissonnat, Devillers, Teillaud and Yvinec4]) to generate tetrahedral meshes for the spatial resolution of diffusion and elasticity via P1 finite elements in both the state and adjoint equations. The meshes used in our simulations consist of roughly
$40$
k tetrahedrons and 48 time steps for the 12 months of regeneration time. The time dependence and couplings in the equations are treated by a semi-implicit ansatz, using only the quantities explicitly that are not available at a current time step due to the couplings of the equations. The ODEs are solved on every element separately, yielding a spatially constant approximation of their solution. Due to the comparatively simple structure of the time-dependent equations, a coarse time stepping can be employed with one temporal increment corresponding to one week of the regeneration process.
4.4. Discussion
We consider two basic set-ups for our optimizaition.
-
A.
$\Omega$
is a simple, cylindrical domain
$B_{10} \times ({-}15,15)$
with a symmetric soft mechanical load of 0.21 kN applied at the two cylinder caps, and zero stress boundary conditions on the mantle. The boundary conditions for
$a_j$
are
$a_j=+1$
,
$j=1,2$
on the cylinder caps, and zero Neumann on the mantle. In this radially symmetric setting, we obtain a nearly radially symmetric optimal scaffold design (Scaffold Architecture A), see Figure 1a. In particular, we see a scaffold with a higher density in the middle region. This is a reasonable outcome, as regenerated bone grows back at the scaffold ends where it is attached to the remaining intact bone tissue. Therefore, the central scaffold region needs to maintain structural integrity for a longer time by itself. The overall shape is very similar to the results obtained by [Reference Poh, Valainis, Bhattacharya, van Griensven and Dondl32] with a one-dimensional model, which is not surprising as our experiment is essentially one dimensional. -
B. We now add an external fixateur to the mechanical environment. Such fixateurs are commonly used in orthopaedic surgery; they usually consist of titanium plates screwed to one side of the healthy remaining bone tissue. In our simulation, for simplicity, we model this fixateur as an additional rod attached to the cylinder. In this rod, the elastic moduli are set to those of titanium with Young’s modulus of 100 GPa and Poisson’s ratio of 0.31. The fixateur domain is excluded from computation of
$b$
,
$c$
,
$a_j$
, and
$\rho$
. The soft load in this case is 0.3 kN. This introduces an asymmetry due to stress shielding. We therefore obtain a non radially symmetric optimal scaffold design (Scaffold Architecture B), see Figure 1b. We see a considerably different distribution compared to the radially symmetric experiment. A higher density in the central part is favourable for the same reason as in the experiment excluding the fixateur; however, in vicinity of the stiff metal plate, a comparatively low scaffold density is optimal. High porosity in this region of the scaffold is beneficial as it increases the mechanical stimulus due to reduced stability and enhances diffusion of bio-active molecules. Both effects promote bone ingrowth in the region close to the fixateur.
Two optimized scaffold densities for different mechanical environments. Architecture A is the result of an optimization routine with a symmetric, essentially one-dimensional setup (disregarding the fixateur), whereas in architecture B the metal fixateur (in gold) is taken into account in the optimization routine in a three-dimensional model.
To illustrate the benefits with respect to stress shielding of the scaffold architecture B over scaffold architecture A, we compare both architectures in a numerical experiment simulating bone regeneration, including external fixation – as an external fixation is virtually always surgically necessary in real-life medical application. Figure 2 shows the strain magnitude distributions at the initial time-point when no bone has regenerated yet. We clearly observe that architecture B mitigates stress shielding in the vicinity of the fixateur in comparison to architecture A. This trend is sustained two months into the regeneration process, as can be observed in Figure 3. We remark that the reduction of stress-shielding is not directly part of the objective function with respect to which the optimization is carried out. Rather, this effect is an implicit favourable consequence of the objective function (31) that advocates for its usage in scaffold design optimization.
This above-mentioned improvement in strain distribution seen for architecture B does indeed lead to an improvement in bone regeneration. While it is difficult to directly compare the objective function as the fixateur dominates the elastic energy, one can see in Figure 4 that the scaffold optimized for the environment with fixateur does show a somewhat improved bone matrix regeneration profile. As expected, the difference is most pronounced in the vicinity of the fixateur.
Comparison of strain magnitudes before bone regeneration. The strain magnitude distributions for the scaffold architectures displayed in Figure 1 are compared in a vertical cut plane before the healing process, i.e., when no bone has regenerated yet. The left picture corresponds to architecture A and the right figure to architecture B of Figure 1.
Comparison of strain magnitudes after two months of bone regeneration. The strain magnitude distributions for the scaffold architectures displayed in Figure 1 are compared in a vertical cut plane two months into the healing process. The left picture corresponds to architecture A and the right figure to architecture B of Figure 1.
Comparison of regenerated bone matrix. The figures show the regenerated bone density after 12 months of regeneration time. The left picture corresponds to architecture A and the right figure to architecture B of Figure 1. One can see that the fixateur-adapted optimized scaffold design leads to improved bone regeneration in the vicinity of the fixateur.
5. Conclusions and future research
We analyzed a three-dimensional, homogenized model for bone growth in the presence of a porous, bio-resorbable scaffold and considered the associated problem of optimal scaffold design. This leads to a PDE-constrained optimization problem for which we proved the existence of an optimal control, i.e., an optimal scaffold density distribution. We presented proof-of-concept numerical experiments illustrating the benefits of a three-dimensional optimization routine. For future work, we propose to use the computational model in detailed numerical simulations and to study the optimized scaffold architectures in vivo.
Funding statement
The authors gratefully acknowledge support from Bundesministerium für Bildung und Forschung (BMBF) within the e:Med programme in the SyMBoD consortium (grant number 01ZX1910C).
Competing interests
The authors declare no competing interests.
Appendix A. proofs of the main results
Lemma 16 (
$C^0(I,H^1(\Omega))$
bound for
$u$
). Let
$\mathbb{C} \in C^0(I,L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s)))$
be uniformly elliptic with ellipticity constant
$\lfloor \mathbb{C} \rfloor$
independent of
$t\in I$
and
$x\in \Omega$
, i.e., it holds
\begin{equation*} \mathbb{C}(t,x)M\,:\,M \geq \lfloor \mathbb{C} \rfloor \lvert M \rvert ^2, \quad \text{for all }M\in \mathcal{L}(\mathcal{M}_s)\text{ and }(t,x)\in I\times \Omega . \end{equation*}
Furthermore, let
$f \in C^0(I,H^1_D(\Omega)^*)$
be a fixed right-hand side. Then the unique solution
$u\in L^2(I,H^1_D(\Omega))$
to
\begin{equation} \iint \mathbb{C}\varepsilon (u)\,:\,\varepsilon ({\cdot})\mathrm dx \mathrm dt = \int _I \langle f,\cdot \rangle _{H_D^{1}(\Omega)}\mathrm dt \quad \text{in } L^2\big(I,H^1_D(\Omega)\big)^* \end{equation}
is in the space
$C^0(I,H^1_D(\Omega))$
and satisfies
\begin{equation*} \lVert u \rVert _{C^0(I,H^1(\Omega))} \leq C\left (\lfloor \mathbb{C}\rfloor , C_{\text{Korn}} \right)\cdot \lVert f \rVert _{C^0(I,H^1_D(\Omega)^*)}. \end{equation*}
Proof. The equation (A1) implies that
$u$
satisfies almost everywhere in
$I$
\begin{equation*} \underbrace {\int _\Omega \mathbb{C}(t)\varepsilon (u(t))\,:\,\varepsilon ({\cdot})\mathrm dx}_{=: \mathcal{T}_tu(t)} = f(t) \quad \text{in }H^1_D(\Omega)^* \end{equation*}
upon applying the isometry
$L^2(I,H^1_D(\Omega))^* \to L^2(I,H^1_D(\Omega)^*)$
to both sides of the equation. Clearly, testing with
$u(t)$
yields, inferring Korn’s inequality,
\begin{equation*} \lfloor \mathbb{C} \rfloor C_{\text{Korn}} \cdot \lVert u \rVert ^2_{H^1_D(\Omega)} \leq \lfloor \mathbb{C} \rfloor \cdot \lVert \varepsilon (u) \rVert _{L^2(\Omega)}^2 \leq \lVert f(t) \rVert _{H^1_D(\Omega)^*}\lVert u(t) \rVert _{H^1_D(\Omega)}. \end{equation*}
Hence,
\begin{equation*} \lVert u(t) \rVert _{H^1_D(\Omega)} \leq \left (\lfloor \mathbb{C} \rfloor C_{\text{Korn}} \right)^{-1} \lVert f(t) \rVert _{H^1_D(\Omega)^*} \leq \left (\lfloor \mathbb{C} \rfloor C_{\text{Korn}} \right)^{-1} \lVert f \rVert _{C^0(I,H^1_D(\Omega)^*)}, \end{equation*}
meaning that the
$H^1(\Omega)$
bound on
$u(t)$
is independent of
$t\in I$
. To show that
$u$
is continuous in time, we compute for
$t,s \in I$
\begin{equation*} f(t) - f(s) = \mathcal{T}_tu(t) - \mathcal{T}_su(s) = \mathcal{T}_t(u(t)-u(s)) + \mathcal{T}_t(u(s)) - \mathcal{T}_su(s). \end{equation*}
Using the coercivity of
$\mathcal{T}_t$
we find
\begin{equation*} \lVert u(t) - u(s) \rVert _{H^1(\Omega)} \leq \frac {1}{\lfloor \mathbb{C} \rfloor C_{\text{Korn}}}\left [ \lVert f(t) - f(s) \rVert _{H^1(\Omega)^*} + \lVert \mathcal{T}_tu(s) - \mathcal{T}_su(s) \rVert _{H^1_D(\Omega)^*} \right ] \end{equation*}
By the assumption
$f\in C^0(I,H^1_D(\Omega)^*)$
it is clear that the first term above tends to zero when
$|t-s|\to 0$
. It remains to estimate
\begin{align*} \lVert \mathcal{T}_tu(s) - \mathcal{T}_su(s) \rVert _{H^1_D(\Omega)^*} & \leq \sup _{\lVert \varphi \rVert _{H^1_D(\Omega)}\leq 1}\int _\Omega \left [ \mathbb{C}(t) - \mathbb{C}(s) \right ]\varepsilon (u(s))\,:\,\varepsilon (\varphi)\mathrm dx \\ & \leq \lVert \mathbb{C}(t) - \mathbb{C}(s) \rVert _{L^\infty (\Omega , \mathcal{M}_s)}\lVert u(s) \rVert _{H^1_D(\Omega)}. \end{align*}
The time-independent bound on
$\lVert u(s) \rVert _{H^1_D(\Omega)}$
and the continuity assumption on
$\mathbb{C}$
imply the assertion.
Lemma 17 (Equi-continuity). Let
$(\rho _k)\subset P$
and
$(b_k)\subset C^0(I, C^0(\Omega))$
be sequences such that
$(b_k, \rho _k)\in \mathrm{dom}(\mathbb C)$
for all
$k\in \mathbb N$
. Furthermore assume that
$(b_k)$
is an equi-continuous sequence in
$C^0(I,C^0(\Omega))$
and
$(f_k)$
is a equi-continuous and bounded sequence in
$C^0(I,H^1_D(\Omega)^*)$
. Assume that
$\mathbb{C}(\rho _k,\sigma , b_k)$
satisfies the assumption 2.2, i.e., in particular, it holds
\begin{equation} \lVert \mathbb{C}(\rho _k,\sigma , b_k)(t) - \mathbb{C}(\rho _k, \sigma , b_k)(s) \rVert _{L^\infty (\Omega , \mathcal{L}(\mathcal{M}_s))} \leq C \lVert b_k(t) - b_k(s) \rVert _{C^0(\Omega)} \end{equation}
for a constant
$C$
that does not depend on the data
$\rho _k\in P$
and
$b_k\in W_{\rho _k}$
and
$t\in I$
. Denote by
$u_k$
the unique solution of
\begin{equation*} \iint \mathbb{C}(\rho _k,\sigma ,b_k)\varepsilon (u_k)\,:\,\varepsilon ({\cdot})\mathrm dx \mathrm dt = \int _I \langle f_k,\cdot \rangle _{H^1_D(\Omega)}\mathrm dt \quad \text{in }L^2\big(I,H^1_D(\Omega)\big)^*. \end{equation*}
Then,
$(u_k)$
lies in
$C^0(I,H^1_D(\Omega))$
and is equi-continuous in this space.
Proof. As the assumptions of Lemma 16 are fulfilled, we know that
$u_k$
is in the space
$C^0(I,H^1_D(\Omega))$
and we need only to establish the equi-continuity. To this end, repeating the equations in Lemma 16 for
$u_k$
instead of
$u$
we arrive at
\begin{align*} \lVert u_k(t) - u_k(s) \rVert _{H^1(\Omega)} &\leq \frac {1}{\lfloor \mathbb{C}_k \rfloor C_{\text{Korn}}}\big[ \lVert f_k(t) - f_k(s) \rVert _{H^1(\Omega)^*} + \lVert \mathbb{C}_k(t) - \mathbb{C}_k(s) \rVert _{L^\infty (\Omega ,\mathcal{L}(\mathcal{M}_s))} \lVert u_k(s) \rVert _{H^1(\Omega)} \big] \\[3pt]&\leq \frac {1}{\lfloor \mathbb{C}_k \rfloor C_{\text{Korn}}}\big[ \lVert f_k(t) - f_k(s) \rVert _{H^1(\Omega)^*} + C\lVert b_k(t) - b_k(s) \rVert _{C^0(\Omega)} \big] \end{align*}
as
$\lVert u_k(t) \rVert _{H^1(\Omega)}$
is bounded uniformly in
$k\in \mathbb{N}$
and
$s\in I$
by Lemma 16 through the boundedness we assumed for
$(f_k)$
. Then, we infer the equi-continuity of
$(f_k)$
and
$(b_k)$
to derive it for
$(u_k)$
.
The following lemma summarizes the main result of [Reference Haller-Dintelmann, Meinlschmidt and Wollner20]. We restrict ourselves to the generality necessary for our application, which, however, is not the most general situation. We refer the reader to [Reference Haller-Dintelmann, Meinlschmidt and Wollner20] for a relaxation concerning boundary regularity, regularity of coefficients and the differential operator.
Lemma 18 (Higher regularity for elliptic systems). Let
$\mathbb{C} \in L^\infty (\Omega , \mathcal{L}(\mathcal{M}_s))$
be uniformly elliptic, i.e., there exists
$\lfloor \mathbb{C} \rfloor \gt 0$
such that
\begin{equation*} \mathbb{C}M\,:\,M \geq \lfloor \mathbb{C} \rfloor |M|^2, \quad \text{for all }M\in \mathcal{M}_s. \end{equation*}
Assume that
$\mathbb{C}_{ijkl}\in C^\alpha (\Omega)$
for a fixed but arbitrary small
$\alpha \gt 0$
. Then, there exists
$\theta = \theta (\alpha) \gt 0$
such that for every
$f\in H^{1-\theta }(\Omega)^*$
the solution
$u\in H^1_D(\Omega)$
to
\begin{equation*} \int _\Omega \mathbb{C}\varepsilon (u)\,:\,\varepsilon ({\cdot})\mathrm dx = f \quad \text{in }H^1_D(\Omega)^* \end{equation*}
is in fact in the space
$H^{1+\theta }(\Omega)$
and we can estimate
\begin{equation*} \lVert u \rVert _{H^{1+\theta }(\Omega)} \leq C\lVert \mathbb{C} \rVert _{C^\alpha (\Omega)}\lVert f \rVert _{H^{1 - \theta }_D(\Omega)^*}, \end{equation*}
where
$C$
does not depend on the concrete form of
$\mathbb C$
.
Proof. This follows from Theorem 1 and Lemma 1 in [Reference Haller-Dintelmann, Meinlschmidt and Wollner20].
The last result we need to establish the relative compactness of
$(u_k)$
in
$C^0(I,H^1_D(\Omega))$
is – little surprisingly – a vector valued version of the Arzelà-Ascoli Theorem which we recall here for convenience.
Theorem 19 (Characterization of relative compactness in
$C^0(K,X)$
spaces). Let
$X$
be a Banach space and
$K$
a compact metric space. Then a set
$\mathcal{F}\subset C^0(K,X)$
is relatively compact if and only if the following two conditions hold:
-
(i) The set
$\mathcal{F}$
is equi-continuous, that is, for all
$t\in K$
and all
$\varepsilon \gt 0$
there exists a neighbourhood
$U(t)\subset K$
such that
\begin{equation*} \sup _{u\in \mathcal{F}}\lVert u(t) - u(s) \rVert _{X} \leq \varepsilon \quad \text{for all }s\in U(t). \end{equation*}
-
(ii) For all
$t\in K$
the set
is relatively compact.
\begin{equation*} \{ u(t)\mid u\in \mathcal{F} \} \subset X \end{equation*}
The focus of the next lemma lies on the a priori estimates for linear parabolic equations.
Lemma 20 (A priori estimate for parabolic evolution equations). Let
$(i,X,H)$
be a Gelfand triple,
$M\,:\,X\to X^*$
a linear coercive operator with coercivity constant
$\lfloor M \rfloor$
, i.e., it holds
\begin{equation*} \langle M a, a\rangle _X \geq \lfloor M \rfloor \lVert a \rVert _X^2, \quad \text{for all }a\in X. \end{equation*}
Let
$I = [0,T]$
denote a time interval and
$f\in L^2(I,X^*)$
a fixed right-hand side. Then there exists a unique solution
$a\in H^1(I,X,X^*)$
to
\begin{equation*} \int _I \langle d_ta, \cdot \rangle \mathrm dt + \int _I\langle Ma, \cdot \rangle _X\mathrm dt = \int _I \langle f,\cdot \rangle _X \mathrm dt, \quad \text{in } L^2(I,X)^* \end{equation*}
Furthermore, the norm of the solution
$a$
can be estimated by
\begin{equation} \lVert a \rVert _{H^1(I,X,X^*)} \leq C\left (\lVert M \rVert _{\mathcal{L}(X,X^*)}, \lfloor M \rfloor ^{-1} \right)\cdot \left (\lVert a(0) \rVert _{H} + \lVert f \rVert _{L^2(I,X^*)}\right)\!, \end{equation}
with
$C$
being monotonously increasing in
$\lVert M \rVert _{\mathcal{L}(X,X^*)}$
and
$\lfloor M \rfloor ^{-1}$
.
Proof. We establish only the estimate (A3); the existence of a solution is the well known maximal regularity result of J. L. Lions, see, for instance, [Reference Ern and Guermond15, Part II, Section 6]. To derive the estimate, we note that by the natural isometry
$L^2(I,X^*) = L^2(I,X)^*$
the function
$a$
satisfies a pointwise almost-everywhere equation in
$X^*$
, namely
\begin{equation*} d_ta(s) + M(a(s)) = f(a(s)) \end{equation*}
which, at time
$s\in I$
, we can test with
$a(s) \in X$
and integrate from
$0$
to
$t$
. Then, we apply the partial integration formula for Gelfand triples and estimate using the coercivity of
$M$
and Young’s inequality
\begin{align*} \frac 12 \lVert a(t) \rVert _H^2 + \lfloor M \rfloor \int _0^t\lVert a(s) \rVert ^2_X\mathrm ds &\leq \frac 12 \lVert a(0) \rVert _H^2 + \lVert f \rVert _{L^2(I,X^*)}\lVert a \rVert _{L^2(I,X)} \\&\leq \frac 12 \lVert a(0) \rVert _H^2 + \frac {1}{2\lfloor M \rfloor }\lVert f \rVert ^2_{L^2([0,t],X^*)} + \frac {\lfloor M \rfloor }{2} \lVert a \rVert ^2_{L^2([0,t],X)}, \end{align*}
which leads to
\begin{equation*} \frac 12 \lVert a(t) \rVert _H^2 + \frac {\lfloor M \rfloor }{2} \int _0^t\lVert a(s) \rVert ^2_X\mathrm ds \leq \frac 12 \lVert a(0) \rVert _H^2 + \frac {1}{2\lfloor M \rfloor }\lVert f \rVert ^2_{L^2([0,t],X^*)}. \end{equation*}
We get by estimating the terms of the left-hand side separately and taking the supremum over
$t\in I$
both
\begin{gather*} \lVert a \rVert _{C^0(I,L^2(\Omega))}^2 \leq \lVert a(0) \rVert _H^2 + \frac {1}{\lfloor M \rfloor }\lVert f \rVert ^2_{L^2(I,X^*)} \quad \text{and} \quad \lVert a \rVert _{L^2(I,X)} \leq \lVert a(0) \rVert _H^2 + \frac {1}{\lfloor M \rfloor ^2}\lVert f \rVert ^2_{L^2(I,X^*)}. \end{gather*}
To estimate the
$L^2(I,X)^*$
norm of
$d_ta$
, we use that
$a$
is the solution of the parabolic equation to estimate
\begin{align*} \lVert d_ta \rVert _{L^2(I,X)^*} &= \sup _{\lVert \varphi \rVert _{L^2(I,X)\leq 1}} \int _I \langle d_ta, \varphi \rangle _X\mathrm dt \\ & \leq \sup _{\lVert \varphi \rVert _{L^2(I,X)\leq 1}} \left [ \int _I |\langle Ma,\varphi \rangle _X|\mathrm dt + \int _I|\langle f,\varphi \rangle _X| \mathrm dt \right ] \\ & \leq \lVert M \rVert _{\mathcal{L(X,X^*)}}\lVert a \rVert _{L^2(I,X)} + \lVert f \rVert _{L^2(I,X^*)}. \end{align*}
If we infer the previous estimates for
$a$
in
$L^2(I,X)$
norm, we can bound
$d_ta$
in
$L^2(I,X)^*$
norm. Combining the considerations for
$a$
and
$d_ta$
lets us bound the
$H^1(I,X,X^*)$
as desired.
Lemma 21 (
$L^p(I,C^\alpha (\Omega))$
bound for
$(a_k)$
). Assume
$\Omega \subset \mathbb{R}^d$
with
$d=1,2,3$
and
$\partial \Omega = \Gamma _N\cup \Gamma _D$
where
$\Omega \cup \Gamma _N$
is Gröger regular. Let
$D\in L^\infty (\Omega ,\mathcal{M}_s)$
be uniformly elliptic with ellipticity constant
$\lfloor D \rfloor \gt 0$
,
$k \gt 0$
,
$f\in L^p(I,L^2(\Omega))$
for a fixed
$p \gt 2$
and
$a_0 \in L^\infty (\Omega)$
some essentially bounded initial condition. Then there exists
$\alpha = \alpha (p)\in (0,1)$
independent of
$D$
and
$f$
such that the solution
$a \in H^1(I,H^1_D(\Omega),$
$H^1_D(\Omega)^*)$
to
\begin{align*} \int _I \langle d_ta,\cdot \rangle _{H^1_D(\Omega)}\mathrm dt + \iint D\nabla a\nabla \cdot + k a\cdot \mathrm dx \mathrm dt &= \iint f\cdot \mathrm dx\mathrm dt \\ a(0) = a_0 \end{align*}
is in
$L^p(I,C^\alpha (\Omega))$
and satisfies the estimate
\begin{equation*} \lVert a \rVert _{L^p(I,C^\alpha (\Omega))} \leq C\left (p,\lfloor D \rfloor , \lVert D \rVert _{L^\infty (\Omega ,\mathcal{M}_s)}\right) \lVert f \rVert _{L^p(I,L^2(\Omega))}. \end{equation*}
Proof. The proof of this Lemma exceeds the scope of this manuscript and is the main result of [Reference Dondl and Zeinhofer14].
Remark 22. We comment on some of the aspects leading to the complexity of the proof of Lemma 21 .
-
(i) The mixed boundary conditions, rough coefficients and the jump initial condition prevent the standard theory from being applicable. If it wasn’t for this roughness, an
$L^2(I,H^2(\Omega))$
result could be derived by standard theory, see, for instance, [Reference Evans16]. -
(ii) Even invoking the theory of abstract parabolic equations as described in [Reference Amann1] does only almost suffice. In fact, combining the results in [Reference Amann1] with [Reference Haller-Dintelmann, Meyer, Rehberg and Schiela21] yields
$L^2(I,C^\alpha (\Omega))$
regularity only if
$a_0$
lies in a suitable trace space for initial conditions. The trace space in this case is
$H^1_D(\Omega)$
and not
$L^\infty (\Omega)$
. -
(iii) The strategy to prove Lemma 21 is therefore to treat the cases
$f = 0$
,
$a(0) = a_0$
and
$f = f$
,
$a(0) = 0$
separately and then to use the superposition principle for linear equations. For details, we refer to [Reference Dondl and Zeinhofer14].
We treat now the cell ODE. We need to establish that a solution in
$W^{1,2}(I,C^\alpha (\Omega))$
exists and is suitably bounded in the data. We already have access to the fact that a long-time solution in
$W^{1,2}(I,C^0(\Omega))$
exists; hence the crucial part is to control the Hölder norm of this solution. This can be done by accessing the solution
$c$
through its formulation as a fixed-point and then estimating its Hölder norm if suitable regularity for the data is given.
Lemma 23.
Let
$a_1$
and
$a_2$
be functions in
$L^4(I,C^\alpha (\Omega))$
with
$a_1, a_2 \geq 0$
. Assume that
$\rho \in C^\alpha (\Omega)$
satisfies
$0 \lt \rho \lt 1$
and
$k_6$
and
$k_7$
are positive constants. Then there exists a solution
$c \in W^{1,2}(I,C^0(\Omega))$
to the equation
\begin{equation*} d_tc = k_6a_1a_2(1+k_7c)\left (1 - \frac {c}{1-\rho } \right)\!, \quad c(0) = 0 \end{equation*}
with
$0\leq c \leq 1$
. Furthermore, we can control the
$\alpha$
-Hölder seminorm of
$c$
in the following way
\begin{equation*} \lfloor c(t) \rfloor _\alpha \leq C\big(\lVert a_1 \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert a_2 \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert \rho \rVert _{C^\alpha (\Omega)} \big), \end{equation*}
with the constant
$C$
being monotone in its arguments.
Proof. The existence of a solution in the space
$W^{1,2}(I,C^0(\Omega))$
was already established in Theorem 3.2 in [Reference Dondl, Poh and Zeinhofer13]. We are only concerned with the control over the Hölder seminorm. To simplify notation, we prove the statement for an ODE of the form
\begin{equation} d_tc = m(1 + c)(1 - \theta c), \quad c(0) = 0 \end{equation}
with
$m \in L^2(I,C^\alpha (\Omega))$
and
$\theta \in C^\alpha (\Omega)$
with
$0 \lt \theta ^{-1}(x) \lt 1$
, which implies that the solution
$c$
to (A4) takes values in the unit interval, i.e.,
$c(t,x)\in [0,1]$
, see Lemma B.5. The existence of
$c$
in
$W^{1,2}(I,C^0(\Omega))$
solving (A4) implies upon applying integrating that
$c(t)$
is given by
\begin{equation*} c(t) = \int _0^t m(s)(1+c(s))(1-\theta c(s))\mathrm ds, \end{equation*}
with the integral being a
$C^0(\Omega)$
valued Bochner integral. As point evaluation at
$x\in \overline {\Omega }$
is continuous and linear from
$C^0(\Omega)$
to
$\mathbb{R}$
, it also holds
\begin{equation*} c(t,x) = \int _0^t m(s,x)(1+c(s,x))(1-\theta (x) c(s,x))\mathrm ds. \end{equation*}
We use the above formula and the triangle inequality to estimate
\begin{align*} |c(t,x) - c(t,y)| &\leq \underbrace {\int _0^t|m(s,x)-m(s,y)|\mathrm ds}_{=: A} + \underbrace {\int _0^t |m(s,x)c(s,x) - m(s,y)c(s,y)| \mathrm ds}_{=: B} \\&+ \underbrace {\int _0^t |m(s,y)\theta (y)c(s,y) - m(s,x)\theta (x)c(s,x)| \mathrm ds}_{=: C} \\&+ \underbrace {\int _0^t |m(s,y)\theta (y)c(s,y)^2 - m(s,x)\theta (x)c(s,x)^2| \mathrm ds}_{=: D}. \end{align*}
For brevity, we set
$\tilde m(t,x) = m(t,x)\theta (x)$
. Inferring that
$c$
takes values in
$[0,1]$
, we claim that the above estimate leads to
\begin{equation} |c(t,x) - c(t,y)| \leq \int _0^t \left (2\lfloor m(s) \rfloor _\alpha + \lVert m(s) \rVert _{C^0(\Omega)}\lfloor c(s) \rfloor _\alpha + 3\lVert \tilde {m}(s) \rVert _{C^0(\Omega)} \lfloor c(s) \rfloor _\alpha + 2\lfloor \tilde {m}(s) \rfloor _\alpha \right)|x-y|^\alpha \mathrm ds. \end{equation}
Dividing by
$|x-y|^\alpha$
and taking the supremum over pairs
$(x,y)\in \overline {\Omega }^2$
with
$x\neq y$
we get
\begin{equation*} \lfloor c(t) \rfloor _\alpha \leq \int _0^t \underbrace {2 \left (\lfloor m(s) \rfloor _\alpha + \lfloor \tilde m(s) \rfloor _\alpha \right)}_{=: \alpha (s)} + \underbrace {\left (\lVert m(s) \rVert _{C^0(\Omega)} + 3\lVert \tilde m(s) \rVert _{C^0(\Omega)} \right)}_{=: \beta (s)} \lfloor c(s) \rfloor _\alpha \mathrm ds. \end{equation*}
Hence, by Grönwall’s lemma we get
\begin{equation*} \lfloor c(t) \rfloor _\alpha \leq \left (1 + \lVert \beta \rVert _{L^1(I)}\exp \left (\lVert \beta \rVert _{L^1(I)}\right) \right)\lVert \alpha \rVert _{L^1(I)} \end{equation*}
with
\begin{align*} \lVert \alpha \rVert _{L^1(I)} &\leq 2 \lVert m \rVert _{L^1(I,C^\alpha (\Omega))} + 2 \lVert \tilde m \rVert _{L^1(I,C^\alpha (\Omega))} \\ \lVert \beta \rVert _{L^1(I)} &\leq \lVert m \rVert _{L^1(I,C^0(\Omega))} + 3 \lVert \tilde m \rVert _{L^1(I,C^0(\Omega))}. \end{align*}
As for bounded intervals, the
$L^2$
norm dominates the
$L^1$
norm, so we are done, given that we still provide the details of the computations that led to (A5). To this end, note that we may estimate
$(A)$
by
\begin{equation*} \int _0^t |m(s,x)-m(s,y)|\mathrm ds \leq \int _0^t \lfloor m(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds. \end{equation*}
Using the triangle inequality and the pointwise properties of
$c$
, we estimate
$(B)$
by
\begin{align*} \int _0^t |m(s,x)c(s,x) - m(s,y)c(s,y)| \mathrm ds &\leq \int _0^t |m(s,x)|\lfloor c(s) \rfloor _\alpha |x - y|^\alpha \mathrm ds + \int _0^t |c(s,y)|\lfloor m(s)\rfloor _{\alpha }|x-y|^\alpha \mathrm ds \\[3pt]&\leq \int _0^t \lVert m(s)\rVert _{C^0(\Omega)} \lfloor c(s) \rfloor _\alpha |x - y|^\alpha \mathrm ds + \int _0^t \lfloor m(s)\rfloor _{\alpha }|x-y|^\alpha \mathrm ds. \end{align*}
Using again the abbreviation
$\tilde m = m\theta$
and noting that
$\tilde m$
has the same regularity as
$m$
, we can estimate the term
$(C)$
in analogy to term
$(B)$
by
\begin{equation*} \int _0^t |m(s,y)\theta (y)c(s,y) - m(s,x)\theta (x)c(s,x)| \mathrm ds \leq \int _0^t \lVert \tilde m(s) \rVert _{C^0(\Omega)}\lfloor c(s) \rfloor _\alpha |x-y|^\alpha + \int _0^t \lfloor \tilde m(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds. \end{equation*}
To estimate
$(D)$
we need to split the term
\begin{equation*} (D) = \underbrace {\int _0^t \tilde |m(s,y)|\,\left | c(s,y)^2 - c(s,x)^2 \right |\mathrm ds}_{=: D_1} + \underbrace {\int _0^t |c(s,x)|^2\left | \tilde m(s,y) - \tilde m(s,x) \right | \mathrm ds}_{=: D_2}. \end{equation*}
Using
$c(s,y)^2 - c(s,x)^2 = c(s,y)(c(s,y)-c(s,x)) + c(s,x)(c(s,y)-c(s,x))$
and
$c(s,x)\in [0,1]$
, we estimate
$(D_1)$
\begin{align*} (D_1) &\leq \int _0^t |\tilde m(s,y)|\, |c(s,y)|\,\lfloor c(s)\rfloor _\alpha |x-y|^\alpha + |\tilde m(s,y)|\,|c(s,x)|\,\lfloor c(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds \\&\leq \int _0^t 2|\tilde m(s,y)|\,\lfloor c(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds \\&\leq \int _0^t 2\lVert \tilde m(s)\rVert _{C^0(\Omega)}\,\lfloor c(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds \end{align*}
and for
$(D_2)$
\begin{equation*} (D_2) \leq \int _0^t |c(s,x)|^2\lfloor \tilde m(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds \leq \int _0^t \lfloor \tilde m(s) \rfloor _\alpha |x-y|^\alpha \mathrm ds. \end{equation*}
Collecting all estimates yields the claim and the proof is complete.
Lemma 24.
Let
$a_1$
and
$a_2$
be functions in
$L^4(I,C^\alpha (\Omega))$
with
$a_1, a_2 \geq 0$
. Assume that
$\rho \in C^\alpha (\Omega)$
satisfies
$0 \lt \rho \lt 1$
and
$k_6$
and
$k_7$
are positive constants. Then there exists a unique solution
$c \in W^{1,2}(I,C^\alpha (\Omega))$
to the equation
\begin{equation*} d_tc = k_6a_1a_2(1+k_7c)\left (1 - \frac {c}{1-\rho } \right)\!, \quad c(0) = 0 \end{equation*}
with
$0\leq c \leq 1$
. Furthermore, we can control the full
$\alpha$
-Hölder norm of
$c$
in the following way
\begin{equation*} \lVert c \rVert _{W^{1,2}(I,C^\alpha (\Omega))} \leq C\left (\lVert a_1 \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert a_2 \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert \rho \rVert _{C^\alpha (\Omega)} \right)\!, \end{equation*}
with the constant
$C$
being monotone in its arguments.
Proof. We use again the notation
\begin{equation*} d_tc = m(1+c)(1-\theta c), \quad c(0)=0, \end{equation*}
where
$m\in L^2(I,C^\alpha (\Omega))$
and
$m(t,x) \geq 0$
and
$\theta \in C^\alpha (\Omega)$
. Thus, the inducing function
$F$
in the sense of Theorem B.2 in [Reference Dondl, Poh and Zeinhofer13] is given by
\begin{equation*} F\,:\,I\times C^\alpha (\Omega) \to C^\alpha (\Omega), \quad F(t,c) = m(t)(1 + c)(1 - \theta c). \end{equation*}
To prove the existence of a unique short-time solution in the space
$W^{1,2}(I_\delta , C^\alpha (\Omega))$
, we need
$F$
to be of Carathéodory regularity. Clearly,
$F(\cdot ,c)\,:\,I\to C^\alpha (\Omega)$
is Bochner measurable as
$m$
is. Furthermore,
$F(t,\cdot)\,:\,C^\alpha (\Omega) \to C^\alpha (\Omega)$
is continuous. This is due to the fact that
$C^\alpha (\Omega)$
is a Banach algebra.
To proceed, we need a boundedness and a Lipschitz condition on bounded subsets of
$C^\alpha (\Omega)$
; compare to Theorem B.2 in [Reference Dondl, Poh and Zeinhofer13]. To this end, let
$B\subset C^\alpha (\Omega)$
be a bounded set. For
$c \in B$
, we estimate
\begin{equation*} \lVert F(t,c) \rVert _{C^\alpha (\Omega)} \leq C \lVert m(t) \rVert _{C^\alpha (\Omega)}\lVert 1 + c \rVert _{C^\alpha (\Omega)}\lVert 1-\theta c \rVert _{C^\alpha (\Omega)} \end{equation*}
The term
$\lVert 1 + c \rVert _{C^\alpha (\Omega)}\lVert 1-\theta c \rVert _{C^\alpha (\Omega)}$
can be bounded in terms of the measure of
$\Omega$
, the assumed boundedness of
$B$
and the
$C^\alpha (\Omega)$
norm of
$\theta$
. Hence, there exists a constant
$C = C(\Omega , \lVert \theta \rVert _{C^\alpha (\Omega)}, B)$
such that
\begin{equation*} \lVert F(t,c) \rVert _{C^\alpha (\Omega)} \leq C\left (\Omega , \lVert \theta \rVert _{C^\alpha (\Omega)},B \right) \lVert m(t) \rVert _{C^\alpha (\Omega)} \end{equation*}
and by assumption, the map
$t\mapsto \lVert m(t) \rVert _{C^\alpha (\Omega)}$
is in
$L^2(I)$
. Now, let
$c$
and
$\bar c \in B$
and look at the differences
\begin{align*} \lVert F(t,c) - F(t,\bar c) \rVert _{C^\alpha (\Omega)} &\leq C \lVert m(t) \rVert _{C^\alpha (\Omega)} \lVert (1+c)(1-\theta c) - (1+\bar c)(1 - \theta \bar c) \rVert _{C^\alpha (\Omega)} \\[2pt]&\leq C\lVert m(t) \rVert _{C^\alpha (\Omega)}\big[ \lVert 1+ \theta \rVert _{C^\alpha (\Omega)} \lVert c - \bar c \rVert _{C^\alpha (\Omega)} + \lVert \theta \rVert _{C^\alpha (\Omega)}\lVert c^2 - \bar c^2 \rVert _{C^\alpha (\Omega)} \big]. \end{align*}
We look at the quadratic term separately
\begin{align*} \lVert c^2 - \bar c^2 \rVert _{C^\alpha (\Omega)} \leq C \lVert c - \bar c \rVert _{C^\alpha (\Omega)}\lVert c + \bar c \rVert _{C^\alpha (\Omega)} \leq C\left (B\right)\lVert c-\bar c \rVert _{C^\alpha (\Omega)}. \end{align*}
Hence, there exists a function
$L_B \in L^2(I)$
such that
\begin{equation*} \lVert F(t,c) - F(t,\bar c) \rVert _{C^\alpha (\Omega)} \leq L_B(t)\lVert c - \bar c \rVert _{C^\alpha (\Omega)}. \end{equation*}
Consulting Theorem B.2 in [Reference Dondl, Poh and Zeinhofer13], the estimates above provide the existence of an interval
$[0,\delta ] = I_\delta$
and a unique function
$c\in W^{1,2}(I_\delta , C^\alpha (\Omega))$
solving the ODE.
To show that the solution can be extended to all of
$I=[0,T]$
, we extend the solution
$c$
to the maximal interval
$[0,t^*)$
of existence. For any
$t_0\in [0,t^*)$
, we set
$c_0 = c(t_0)$
and consider the initial value problem
\begin{equation*} d_tc = k_6a_1a_2(1+k_7c)\left (1 - \frac {c}{1-\rho } \right)\!, c(t_0) = c_0. \end{equation*}
Then this has a unique solution in
$W^{1,2}([t_0-\tilde \delta , t_0 + \tilde \delta ]\cap I, C^\alpha (\Omega))$
for some suitable
$\tilde \delta$
. In fact,
$\tilde \delta$
depends on the
$L^2(I,C^\alpha (\Omega))$
norm of
$a_1a_2$
, the
$C^\alpha (\Omega)$
norm of
$(1-\rho)^{-1}$
and the
$L^2(I,C^\alpha (\Omega))$
norm of
$c$
on
$[0,t^*)$
. This implies that
$\tilde \delta$
does not depend on the position of
$t_0\in [0,t^*)$
and thus
$t^*=T$
and the interval
$[0,t^*)$
can be closed.
Finally, the promised bound on the
$W^{1,2}(I,C^\alpha (\Omega))$
norm of
$c$
is easily established using
$c(t,x)\in [0,1]$
and the estimate on the Hölder seminorm of Lemma 23.
We show now how to establish the existence of solutions to the bone ODE in the space
$W^{1,2}(I,C^\alpha (\Omega))$
. Furthermore, we show that the
$W^{1,2}(I,C^\alpha (\Omega))$
norm of such solutions can be bounded, given bounded data in the right spaces. This can in principle be done by the same arguments as for the cell equation; however, the bone ODE is linear and thus we can use more elegant approaches.
Lemma 25.
Let
$X$
be a Banach algebra and denote by
$C_X \gt 0$
the norm of its multiplication and assume that
$p \gt 1$
. By
$W^{1,p}_0(I,X)$
, we denote the vector-valued Sobolev space with vanishing initial conditions. For a function
$m\in L^p(I,X)$
, we define the multiplication operator
\begin{equation*} M\,:\, C^0(I,X) \to L^p(I,X), \quad Mv = t\mapsto m(t)v(t). \end{equation*}
Then the map
\begin{equation*} d_t + M\,:\, W^{1,p}_0(I,X) \to L^p(I,X), \quad v\mapsto d_tv + Mv \end{equation*}
is a linear homeomorphism. Furthermore, given a right-hand side
$f\in L^p(I,X)$
, we may bound the solution
$v$
to
$d_tv + Mv = f$
in the following way
\begin{equation*} \lVert v \rVert _{W^{1,p}_0(I,X)} \leq C \left (|I|, \lVert m \rVert _{L^p(I,X)}, C_X \right)\lVert f \rVert _{L^p(I,X)}, \end{equation*}
i.e., the norm of
$v$
does only depend on
$f$
and
$m$
measured in
$L^p(I,X)$
norm and the constant
$C$
is monotone in these quantities.
Proof. The continuity and linearity of the map
$d_t + M$
is clear. Its bijectivity follows as an application of Theorem B.2 in [Reference Dondl, Poh and Zeinhofer13]. To this end, note that the inducing function
$F\,:\,I\times X \to X$
of Theorem B.2 in [Reference Dondl, Poh and Zeinhofer13] is given by
\begin{equation*} F\,:\,I\times X \to X, \quad F(t,x) = m(t)x. \end{equation*}
This is clearly a Carathéodory function and it holds for
$x,y \in X$
\begin{equation*} \lVert F(t,x) - F(t,y) \rVert _X \leq C \lVert m(t) \rVert _X\lVert x-y \rVert _X. \end{equation*}
The function
$C\lVert m({\cdot}) \rVert _X$
is in
$L^p(I)$
with
$p \gt 1$
and therefore the existence of a unique solution
$v \in W^{1,p}_0(I,X)$
is established. To provide the bound, we employ Grönwall’s inequality. Note that, by the fundamental theorem, the solution
$v$
satisfies the integral identity
\begin{equation*} v(t) = \int _0^t f(s) - m(s)v(s)\mathrm ds \end{equation*}
and consequently the estimate
\begin{equation*} \lVert v(t) \rVert _X \leq \int _0^t \lVert f(s) \rVert _X + C_X\lVert m(s) \rVert _X\lVert v(s) \rVert _X\mathrm ds. \end{equation*}
Using Grönwall’s inequality yields
\begin{equation*} \lVert v(t) \rVert _X \leq [ 1 + C_X\lVert m \rVert _{L^1(I,X)}\exp (C_X \lVert m \rVert _{L^1(I,X)} ) ]\cdot \lVert f \rVert _{L^1(I,X)}. \end{equation*}
Clearly, this implies a bound in
$C^0(I,X)$
norm for
$v$
of the form
\begin{equation*} \lVert v \rVert _{C^0(I,X)}\leq C (\lVert m \rVert _{L^1(I,X)}, C_X )\lVert f \rVert _{L^1(I,X)} \end{equation*}
and consequently also in
$L^p(I,X)$
. To bound
$d_tv$
, we use the equation satisfied by
$v$
and estimate
\begin{align*} \lVert d_tv \rVert _{L^p(I,X)} &= \lVert f - Mv \rVert _{L^p(I,X)} \\ &\leq \lVert f \rVert _{L^p(I,X)} + \lVert v \rVert _{C^0(I,X)}\lVert m \rVert _{L^p(I,X)} \\ &\leq C (|I|, \lVert m \rVert _{L^p(I,X)}, C_X )\lVert f \rVert _{L^p(I,X)}. \end{align*}
Lemma 26.
Assume
$(\rho _k)$
is a minimizing sequence for
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
. Then properties
$(A1)$
–
$(A5)$
hold.
Proof. We begin with
$(A1)$
. The regularizing term
$\eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
leads to a
$H^2(\Omega)$
bound for any minimizing sequence
$(\rho _k)$
of
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
, as
$\,\hat{\!J}$
is bounded from below. Then there exists a subsequence (not relabelled) with
\begin{equation*} \rho _k \rightharpoonup \rho ^*\quad \text{in }H^2(\Omega). \end{equation*}
Using the compactness of the embedding
$H^2(\Omega)\hookrightarrow \hookrightarrow C^0(\Omega)$
, we get
\begin{equation*} \rho _k \to \rho ^* \quad \text{in }C^0(\Omega) \end{equation*}
and inferring the closedness of
$P$
in
$C^0(\Omega)$
yields
$\rho ^* \in P$
as desired.
We provide first a weaker statement than
$(A2)$
. Namely, we prove that
$(u_k)$
is bounded uniformly in
$C^0(I,H^1(\Omega))$
. At the end of the proof, we can show the full validity of
$(A2)$
. In Setting 2.2 we assumed that the boundary conditions satisfied by
$u_k$
are
\begin{gather*} \mathbb{C}(\rho _k,\sigma ,b_k)\varepsilon (u_k)\cdot \eta = g_N, \quad u_k=u_D \end{gather*}
on
$\Gamma _N$
and
$\Gamma _D$
, respectively, where
$g_N \in C^0(I,L^2(\partial \Omega)) \subset C^0(I,H^{1/2}(\partial \Omega)^*)$
and
$u_D\in C^0(I,H^{1+\theta }(\Omega))$
implying that
$(u_D)_{|\Gamma _D}\in H^{1/2}(\partial \Omega)$
. The unique solutions
$\tilde u_k \in L^2(I,H^1_D(\Omega))$
to
\begin{equation} \iint \!\mathbb{C}(\rho _k, \sigma , b_k)\varepsilon (\tilde u_k)\,:\,\varepsilon ({\cdot})\mathrm dx\mathrm dt = \underbrace {\int _I\langle g_N, \cdot \rangle _{H^{1/2}(\partial \Omega)} \mathrm dt - \iint \mathbb{C}(\rho _k, \sigma , b_k)\varepsilon (u_D)\,:\,\varepsilon ({\cdot})\mathrm dx\mathrm dt}_{=: f_k \in L^2\left(I,H^1_D(\Omega)\right)^*} \quad \text{in }L^2\big(I,H^1_D(\Omega)\big)^* \end{equation}
have therefore right-hand sides
$f_k$
that can be interpreted as elements of
$C^0(I,H^1_D(\Omega)^*)$
. This is due to the assumption
$g_N \in C^0(I,H^{1/2}(\partial \Omega)^*)$
and a standard computation that shows that the map
\begin{equation*} t \mapsto \int _\Omega \mathbb{C}(\rho _k,\sigma , b_k)(t)\varepsilon (u_D(t))\,:\,\varepsilon ({\cdot})\mathrm dx \end{equation*}
is in
$C^0(I,H^1_D(\Omega)^*)$
. Hence Lemma 16 is applicable and shows that
\begin{equation*} \lVert \tilde u_k \rVert _{C^0(I,H^1(\Omega))} \leq C\left (\lfloor \mathbb{C}(\rho _k,\sigma ,b_k) \rfloor , C_{\text{Korn}} \right) \lVert f_k \rVert _{C^0(I,H^1_D(\Omega)^*)}. \end{equation*}
As
$\mathbb{C}$
is coercive and essentially bounded, this estimate can be made independent of
$k\in \mathbb{N}$
. Clearly, we have not yet proven
$(A2)$
but will first continue with the other assumptions.
We are concerned with
$(A3)$
now, which is an application of Lemma 20. The corresponding Gelfand triple is
$({\mathrm{Id}},H^1_D(\Omega),L^2(\Omega))$
and the operators
$M_k$
are
\begin{equation*} M_k\,:\,H^1_D(\Omega) \to H^1_D(\Omega)^*, \quad M_ka = \int _\Omega D(\rho _k)\nabla a\nabla \cdot + k_3a\cdot \mathrm dx. \end{equation*}
The coercivity constant of
$M_k$
can be estimated from below independently of
$(\rho _k)$
by
\begin{equation*} \lfloor M_k \rfloor \geq \min (\lfloor D(\rho _k) \rfloor , k_3). \end{equation*}
On the other hand, the operator norm of
$M_k$
can be estimated to
\begin{equation*} \lVert M_k \rVert = \sup _{\lVert a \rVert \leq 1, \lVert \varphi \rVert \leq 1}\int _\Omega D(\rho _k)\nabla a \nabla \varphi + k_3 a\varphi \mathrm dx \leq \lVert D(\rho _k) \rVert _{L^\infty (\Omega ,\mathcal{M}_s)} + k_3. \end{equation*}
The right-hand sides of the equation are given by
\begin{equation*} f_k = \iint (k_2S(\varepsilon (u_k)) c_k - k_3) \cdot \mathrm dx \mathrm dt, \end{equation*}
consequently their norm can be estimated
\begin{align*} \lVert f_k \rVert _{L^2(I,H^1_D(\Omega)^*)} & \leq \lVert f_k \rVert _{L^2(I,L^2(\Omega)^*)} \\ &\leq \left [ \iint (k_2S(\varepsilon (u_k)) c_k - k_3) \right ]^{1/2} \\ &\leq k_2\lVert c_k \rVert _{C^0(I\times \Omega)}\lVert u_k \rVert _{L^2\left(I,H^1_D(\Omega)\right)} + |I\times \Omega |^{1/2}k_3. \end{align*}
This is uniformly bound in
$k\in \mathbb{N}$
by the bound on
$(u_k)$
and the pointwise properties of
$c_k$
, i.e.,
$0\leq c_k \leq 1$
. Thus we apply Lemma 20 to obtain
\begin{align*} \lVert a_{i, k} \rVert _{H^1(I,H^1_D(\Omega), H^1_D(\Omega)^*)} &\leq C(\lVert M_k \rVert , \lfloor M_k\rfloor ^{-1})\cdot \big(\lVert a_{i, k} (0) \rVert _{L^2(\Omega)} + \lVert f_k \rVert _{L^2\left(I,H^1_D(\Omega)\right)} \big) \\&\leq C\left (\lVert D(\rho _k) \rVert _{L^\infty (\Omega ,\mathcal{M}_s)}, \lfloor D(\rho _k) \rfloor ^{-1}, k_{3,i} \right)\cdot \big(\lVert a_{i, k}(0) \rVert _{L^2(\Omega)} + \lVert f_k \rVert _{L^2\left(I,H^1_D(\Omega)\right)} \big) \end{align*}
By the ellipticity and boundedness of
$D$
, the constant initial conditions
$ a_{i, k}(0)\equiv 1$
and the estimate for
$u_k$
, we see that
$(a_{i,k})$
are bounded uniformly in
$k\in \mathbb{N}$
. Using the reflexivity of the Hilbert space
$H^1(I,H^1_D(\Omega), H^1_D(\Omega)^*)$
to produce a weakly convergent subsequence with limit
$a_i^*$
we proved
$(A3)$
.
We proceed with
$(A4)$
and aim to apply Lemma 25 to the ODE
\begin{equation} d_tb_k = k_4 a_{1, k} \left (1 + \frac {b_k}{1-\rho _k} \right)\!, \quad b_k(0) = 0, \end{equation}
with
$X = C^\alpha (\Omega)$
and
$p = 2$
. To this end we rearrange (A7) to
\begin{equation*} d_tb_k - \frac {k_4 a_{1, k} }{1-\rho _k}b_k = k_4 a_{1, k} , \end{equation*}
thus
\begin{gather*} m_k = \frac {k_4 a_{1, k} }{1 - \rho _k}\quad \text{and}\quad f_k = k_4 a_{1, k} \end{gather*}
in the notation of Lemma 25. Due to the embedding
$H^2(\Omega)\hookrightarrow \hookrightarrow C^\alpha (\Omega)$
in three spatial dimensions, we get
$\rho _k\in C^\alpha (\Omega)$
and also
$(1 - \rho _k)^{-1}\in C^\alpha (\Omega)$
with a uniform bound in Hölder norm. Then, Lemma 25 guarantees that
$(a_{1, k})$
is bounded uniformly in
$L^2(I,C^\alpha (\Omega))$
which implies such a bound for
$(m_k)$
and
$(f_k)$
. We may therefore use Lemma 25 to obtain
\begin{equation*} \lVert b_k \rVert _{W^{1,2}(I,C^\alpha (\Omega))} \leq C\left (I, C_{C^\alpha (\Omega)}, \lVert m_k \rVert _{L^2(I,C^\alpha (\Omega))} \right)\lVert f_k \rVert _{L^2(I,C^\alpha (\Omega))} \end{equation*}
and guarantee that the bound is independent of
$k\in \mathbb{N}$
. Furthermore, we have the compact embedding
\begin{equation*} W^{1,2}(I,C^\alpha (\Omega))\hookrightarrow \hookrightarrow C^0(I\times \Omega). \end{equation*}
This yields the relative compactness of
$(b_k)$
in
$C^0(I\times \Omega)$
and thus the existence of
$b^*\in C^0(I\times \Omega)$
and a (not relabeled) subsequence with
\begin{equation*} b_k \to b^*. \end{equation*}
To provide the existence of
$c^*$
and a subsequence
$c_k\to c^*$
in
$C^0(I\times \Omega)$
, we note that Lemma 24 provides a bound of the
$W^{1,2}(I,C^\alpha (\Omega))$
norm of
$c_k$
of the form
\begin{equation*} \lVert c_k \rVert _{W^{1,2}(I,C^\alpha (\Omega))} \leq C \left (\lVert a_{1, k} \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert a_{2, k} \rVert _{L^2(I,C^\alpha (\Omega))}, \lVert \rho _k \rVert _{C^\alpha (\Omega)} \right) \end{equation*}
with
$C$
being increasing in its arguments. As we proved suitable bounds for
$(a_{1, k})$
,
$(a_{2, k})$
and
$(\rho _k)$
this yields a
$W^{1,2}(I, C^\alpha (\Omega))$
bound for
$(c_k)$
that does not depend on
$k\in \mathbb{N}$
. Therefore, the sequence
$(c_k)$
is relatively compact in
$C^0(I\times \Omega)$
and
$(A5)$
follows.
We are still left with showing the existence of a subsequence of
$(u_k)$
and a function
$u^*\in C^0(I,H^1(\Omega))$
such that
\begin{equation*} u_k \to u^*\quad \text{in }C^0(I,H^1(\Omega)). \end{equation*}
To this end, note that we have established that
$(b_k)$
is relatively compact in
$C^0(I,C^0(\Omega))$
. Hence, applying the Arzelà-Ascoli Theorem,
$(b_k)$
is equi-continuous as well. Now, going back to (A6) we can easily compute that the sequence
$(f_k)$
is indeed equi-continuous in
$C^0(I,H^1_D(\Omega)^*)$
. Again, this is essentially due to the equi-continuity of
$(b_k)$
in the space
$C^0(I,C^0(\Omega))$
. Hence, applying Lemma 17 yields the equi-continuity of
$(u_k)$
in
$C^0(I,H^1(\Omega))$
. In view of the Arzelà-Ascoli Theorem, we still need to show that the sets
\begin{equation*} \{ u_k(t) \in H^1(\Omega)\mid k\in \mathbb{N} \} \end{equation*}
are relatively compact in
$H^1(\Omega)$
. This can be established by looking at the equation satisfied by
$u_k(t)$
for every fixed
$t\in I$
and applying Lemma 18. Indeed,
$u_k(t) =\tilde u_k(t) + u_D(t)$
satisfies
\begin{equation*} \int _\Omega \mathbb{C}(\rho _k,\sigma ,b_k)(t)\varepsilon (\tilde u_k(t))\,:\,\varepsilon ({\cdot})\mathrm dx = \int _{\partial \Omega }g_N(t)\cdot \mathrm ds - \int _\Omega \mathbb{C}(\rho _k,\sigma ,b_k)(t)\varepsilon (u_D(t))\,:\,\varepsilon ({\cdot})\mathrm dx. \end{equation*}
The assumptions on
$g_N$
and
$u_D$
guarantee that the right-hand side lies in
$H^{1-\theta }(\Omega)^*$
and the Hölder regularity established for
$(b_k)$
, i.e.,
$(b_k(t))\subset C^\alpha (\Omega)$
allows to deduce the
$H^{1+\theta }(\Omega)$
regularity for
$u_k(t)$
. Additionally, the uniform
$W^{1,2}(I,C^\alpha (\Omega))$
bound for the sequence
$(b_k)$
established before yields a bound for the
$C^\alpha (\Omega)$
norm of
$(b_k(t))$
and also the coefficients of
$\mathbb{C}(\rho _k,\sigma ,b_k)$
via assumption (17). Collecting these bounds in fact implies that
\begin{equation*} \sup _{k\in \mathbb{N}}\lVert u_k(t) \rVert _{H^{1+\theta }(\Omega)} \leq C \end{equation*}
and invoking the compactness result of Rellich, which states that
\begin{equation*} H^{1+\theta }(\Omega)\hookrightarrow \hookrightarrow H^1(\Omega) \end{equation*}
we can conclude the missing piece in order to apply the Arzelà-Ascoli Theorem to
$(u_k)$
in the space
$C^0(I,H^1_D(\Omega))$
. This eventually guarantees the validity of assumption
$(A2)$
.
We can now conclude the section by providing the proofs of Theorem 11 and Corollary 12.
Proof of Theorem
11. Let
$(\rho _k)\subset P$
be a minimizing sequence for
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)}$
. Lemma 26 shows that the assumptions
$(A1){-}(A5)$
hold and Proposition 14 shows that this leads to the existence of an accumulation point
$\rho ^*\in P$
of the sequence
$(\rho _k)$
which is a minimizer of
$\,\hat{\!J} + \eta \lVert \cdot \rVert _{H^2(\Omega)}^2$
.
Proof of Corollary
12. Let
$(\rho _k)\subset H^2(\Omega)$
be a minimizing sequence for
$\,\hat{\!J} + \eta \lVert \cdot \rVert ^2_{H^2(\Omega)} + \mathcal{K}$
. Revisiting the proof of Proposition 14 shows that the additional term
$\mathcal{K}$
does not lead to complications in the lower semicontinuity as it is assumed to be continuous on
$C^0(I\times \Omega)$
, i.e., a compact perturbation. Furthermore, as we assumed that
$\mathcal{K}$
takes non-negative values only, the coercivity of the objective function is also not violated through the addition of
$\mathcal{K}$
.
Open access
