An Efficient Model for Scaffold-Mediated Bone Regeneration
AAN EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION
PATRICK DONDL, PATRINA S. P. POH, AND MARIUS ZEINHOFER
Abstract.
We present a three dimensional, time dependent model for bone regeneration in the pres-ence of porous scaffolds to bridge critical size bone defects. Our approach uses homogenized quanti-ties, thus drastically reducing computational cost compared to models resolving the microstructuralscale of the scaffold. Using abstract functional relationships instead of concrete effective materialproperties, our model can incorporate the homogenized material tensors for a large class of scaf-fold microstructure designs. We prove an existence and uniqueness theorem for solutions based ona fixed point argument. We include the cases of mixed boundary conditions and multiple, interact-ing signalling molecules, both being important for application. Furthermore we present numericalsimulations showing good agreement with experimental findings.
Contents
1. Introduction 21.1. Scaffold Mediated Bone Growth 21.2. The System of Equations 41.3. Concrete Examples. 52. Mathematical Formulation 62.1. The Domain 62.2. Admissible Data 72.3. The Elastic Equation 72.4. Diffusion Equations 82.5. Ordinary Differential Equations 93. Existence and Uniqueness 104. Numerical Experiments 154.1. Model Parameters 164.2. Numerical Implementations 174.3. Discussion of Numerical Simulations 18Acknowledgements 19Appendix A. Properties of Diffusion Equation 19A.1. Regularity Results 19A.2. Boundary Regularity 20A.3. Proofs of the Regularity Results 20Appendix B. Ordinary Differential Equations 23References 27
Date : January 25, 2021.2020
Mathematics Subject Classification.
Key words and phrases.
Mathematical modeling, tissue engineering, scaffold mediated bone regeneration, coupledPDE systems, mixed boundary conditions. a r X i v : . [ m a t h . A P ] J a n PATRICK DONDL, PATRINA S. P. POH, AND MARIUS ZEINHOFER Introduction
In this work, we are concerned with the development and well-posedness of a simple and efficient modelfor bone regeneration in the presence of a bioresorbable porous scaffold. The essential processes arean interplay between the mechanical and biological environment which we model by a coupled systemof PDEs and ODEs. The mechanical environment is represented by a linear elastic equation and thebiological environment through reaction-diffusion equations as well as as logistic ODEs, modelling sig-nalling molecules and cells/bone respectively. Material properties are incorporated using homogenizedquantities not resolving any scaffold microstructure. This makes the model efficient in computations,thus suitable as a forward equation in optimization algorithms and opening up the possibility of patientspecific scaffold design in the sense of precision medicine.We analyze the model mathematically, proving well-posedness. We stress that we allow data that isrealistic for applications, i.e., non-smooth domains and mixed Dirichlet-Neumann boundary conditions.The article is organized as follows. Next, we give an introduction into tissue engineering for the treat-ment severe bone defects and present our computational model. Then, we discuss its weak formulationin section 2 and prove an existence and uniqueness result in section 3. Finally, numerical simulationsare presented in 4. Appendix A is concerned with regularity results for Dirichlet-Neumann boundaryvalue problems and Appendix B contains results on Banach space valued ODEs. We include the latterbecause, even though the results are folklore, we are not aware of any references.1.1.
Scaffold Mediated Bone Growth.
The regeneration and restoration of skeletal functions ofcritical-sized bone defects ( >
25 mm) are very challenging despite a multitude of treatment options[34]. The main problem is the phenomenon of non-union where the bone defect fails to become bridgedafter > . £
320 million annually [47]. Moreover, the risk of non-unionincreases drastically with comorbidities such as diabetes as in this case the regenerative capability ofbone tissue is compromised [32].Critical-sized defects may not heal and require in-depth planning of their treatment. Currently usedtherapeutic approaches include bone grafting, distraction osteogenesis, and the so-called “Masquelet”technique, in which a periosteal membrane is formed to induce bone defect healing [34]. Despite havinga general guideline for treatment of critical-sized bone defects, healing outcomes vary highly, depen-dent on the site and size of the defect and patient-related aspects, e.g., age, lifestyle and comorbidmetabolic/systemic disorders [42].Over the years, research illustrated the potential of using porous, possibly bio-resorbable supportstructures, so-called scaffolds, as supporting devices to promote bone defect regeneration. Initially, ascaffold is placed in the defect site, acting as a temporary support structure allowing for vascular-ization while guiding new bone formation. This has recently shown promising results in vivo and inclinical cases, for example [36] showed that the architecture of the scaffold can guide the endochondralhealing of bone defects in rats. In this study, collagen-based scaffolds with cylindrical pores alignedalong the principle stress axis were used. In [14, 35], 3D-printed scaffolds made from a composite ofpolycaprolactone (PCL, a slowly degrading, bio-resorbable synthetic thermoplastic) and β -tricalciumphosphate ( β -TCP) were used in an ovine experiment. In the studies [36, 14, 35] no relevant bridging ofthe bone defect was achieved without the addition of exogenous growth factors or cells. However, [38]illustrated that clinically relevant bone formation for scaffold mediated bone regeneration is possiblewithout exogenous growth factors. In this experiment a 3D-printed titanium scaffold with optimizedmechanobiological properties was used and displayed clinically relevant functional bridging of a majorbone defect in a large animal model. Concluding, these studies [36, 14, 35, 38] indicate the possibilityof using a scaffold-mediated bone growth approach for critical-size bone defect healing. Furthermorethey indicate that the design and choice of materials are critical questions not yet fully understood.There are several objectives to be considered when designing a scaffold, such as (a) the porosity, poresize and shape, influencing cell proliferation and differentiation as well as the vascularization process; N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 3 (b) the overall stability and elastic properties guaranteeing a proper transfer of loads, as mechanicalstimulus is indispensable for bone growth; (c) patient specific information such as reduced bone healingcapacities, caused for example by diabetes [32]. Therefore, the patient dependent optimal scaffolddesign is of fundamental importance and with the advent of additive manufacturing technologies theproduction of personalized scaffolds is – in theory – fully feasible.However, the design of scaffolds has been dominated by trial-and-error approaches – modifying anexisting scaffold architecture based on experimental outcomes, a very costly workflow unsuitable forpatient specific design. Over the years, with the help of evolving computer aided design tools, topologyoptimization techniques have shown potential to address the optimal design question computationally.This strategy has already been applied to design scaffolds meeting elastic optimality conditions with agiven porosity or fluid permeability [17, 16, 30, 23, 10, 25, 50, 19]. Yet, a common limitation to thesemodels is that they do not resolve the time dependence of the bone regeneration process, as scaffoldmediated bone regeneration crucially depends on the varying elastic moduli over time.Highly accurate, fine scale models for bone formation exist (see, e.g., [28, 43, 2, 11]). A central issue inmost such micro-scale models is that their use in optimization routines for scaffold design is impededby too high computational cost. Ideally, a bone regeneration scaffold design should be patient specific,i.e., depend on the individual patient’s defect site and its biomechanical loading conditions, geometry,and regenerative ability as influenced by, e.g., comorbitities such as type 2 diabetes mellitus. Such anoptimization of course relies on the availability of highly efficient models for bone regeneration thatnevertheless take into account mechanics and biological signalling.Based on a previous, one-dimensional study [39], we thus propose a model based on homogenizedquantities suitable for scaffold optimization in the sense of the first step in the “Shape Optimizationby the Homogenization Method” [3]. This means that our model does not resolve the micro-structureof the scaffold design, but uses coarse-grained values instead. In a scaffold based on a unit cell design,the scaffold volume fraction (or equivalently, the porosity) changes on a larger length-scale than theunit cell design. We use this fact to simplify our model, working with meso-scale averages of thevolume fraction instead of the precise micro-structure. Likewise, the other quantities of the model canbe viewed as locally averaged values. However, it should be made clear that using such an approachimplies that only the averaged quantities can be tracked over the regeneration process and no predictionon how the micro-structure changes over time can be made. Rather, this is required as an input toprovide the correct homogenized material properties. Our central assumption is that one can describethe time-evolution of the homogenized quantities in terms of their averages at the initial time-point.Compared to the aforementioned one-dimensional approach, our model can resolve important issuessuch as bone mass loss due to stress shielding in orthopaedic implants, see section 4 for an explicitexample.As our model is designed for computational efficiency we include only key events in the course of thebone healing process. We keep track of the mechanical environment at every point in time and space,depending on the current state of bone formation and scaffold degradation in terms of its molecularweight. Here we focus on additively manufactured scaffolds made out of PCL, a very promising materialfor this specific application. Of course, extensions to other materials (e.g., non-degrading titanium) arepossible. The biological environment is represented via a concentration of endogenous angiogenic andosteoinductive factors (e.g., intrinsic growth factors/cytokines) which we call bio-active or signallingmolecules and a concentration of osteoblasts, a type of bone forming cell. The coupling of the mechan-ical and biological properties is assumed to be driven through the local strain caused by mechanicalloading of the scaffold-bone composite, i.e., mechanical loading leads to stimulus for the biologicalenvironment which in turn leads to bone growth and hence changes the mechanical properties.This results in a coupled system of evolution equations composed of a linear elastic equilibrium equa-tion for every point in time, diffusion equations for the bio-active molecules and ordinary differentialequations for the concentration of osteoblasts and the volume fraction of bone. As our main mathemat-ical result we prove that this system admits a unique solution in a certain weak sense, see Theorem 3.2.This shows that our model is well-posed, a necessary requirement for a reasonable biological model.
PATRICK DONDL, PATRINA S. P. POH, AND MARIUS ZEINHOFER
The strategy used to prove Theorem 3.2 is to apply a fixed-point theorem on a map associated withthe coupled system of equations, see the beginning of section 3 for a precise description. The maindifficulty we encounter is the notorious low regularity of mixed Dirichlet-Neumann boundary valueproblems [44, 26, 21] which one is forced to consider when one desires to allow for realistic boundaryconditions See Appendix A where we collect results from the literature that are helpful in our case.As our main focus lies on the existence and uniqueness results, we do not use concrete homogenizedtensors in the equations, but abstract functional relationships. This has the advantage of proving theresult for a wide class of imaginable scaffold architectures at once. The concrete micro-structure canthen be taken into account when one performs numerical simulations. In the same spirit we keepthe rest of the equations abstract, preferring functional relationships over concrete formulas. Thisconstitutes also a perspective for future research: derive concrete homogenized quantities for certainscaffold details, compare the outcome to experimental results, and employ the model in an optimizationroutine analogous to the one presented in [39]. The 3-dimensionality of the model makes an optimizationof the scaffold porosity considerably more challenging from a numerical viewpoint – but due to theefficient, homogenized, model it is within reach to provide patient specific optimal scaffold designs thatdepend on the individual’s defect site and geometry, as well as their regeneration capacity.1.2.
The System of Equations.
Let Ω ⊂ R be the domain of computation, i.e., the bone defect site,and let I = [0 , T ] be some finite time interval. On the defect site we keep track of the local scaffoldvolume fraction called ρ ( x ), with x ∈ Ω. Equivalently, the relation to the local scaffold porosity θ is given by θ ( x ) = 1 − ρ ( x ), but we work with ρ exclusively. Note that we do not assume a timedependency for ρ as experimental findings [37] have shown that, in the time-window relevant for us,PCL degrades via bulk erosion. However, the molecular mass decreases and we keep track of this byintroducing the exponential decay σ ( t ) = e − k t , making the product ρ ( x ) · σ ( t ) the quantity encodingthe mechanical properties of PCL over time and space. Furthermore, we denote the local bone densityby b ( t, x ) and the three quantities b, σ and ρ together determine the mechanical material propertiesof the bone-scaffold composite. We model this composite in the linear elastic regime using an elastictensor C ( ρ, σ, b ) to capture the material properties.In the spirit of the homogenization approach we assume little on the concrete properties of this tensor,in particular we do not assume isotropy. For a particular choice of micro-structure C ( ρ, σ, b ) can bemade explicit. In order to quantify the elastic stimulus throughout the bone-scaffold composite weintroduce a displacement field u ( t, x ) satisfying the equation of mechanical equilibrium (1.1). Thecorresponding strain is denoted by ε ( u ), with ε ( u ) = ( Du + D T u ) the symmetrized derivative.For the biological environment we introduce N bio-active molecules denoted by a ( t, x ) , . . . , a N ( t, x ),these are endogenous angiogenic and osteoinductive factors which we assume to diffuse dependingon the scaffold density ρ . This is captured by D i ( ρ ) in the equation (1.2) and is left as an abstractfunctional relationship for the same reasoning as the elastic tensor. Furthermore, we assume the bio-active molecules to decay at a certain rate and to be produced in the presence of strain and a localdensity of specific cells (e.g., osteoblasts) which we denote by c ( t, x ). The essential quantity for theproduction of bio-active molecules is | ε ( u ) | δ , where | · | δ is a functional relationship which we proposeto view as a usual Euclidean norm or a truncated version thereof, see also (2.12). The concentrations ofbio-active molecules are normalized to unity in healthy tissue and the choice of decay and productionrate should reflect this in a concrete simulation.Equation (1.3) governing the production of bone forming cells (here: osteoblasts) is modeled by lo-gistic growth and a functional relationship H ( a , . . . , a N , c, b ) allowing driving factors for osteoblastproduction to be the concentrations of bio-active molecules (causing differentiation of stem cells toosteoblasts), the proliferation of osteoblasts and the maturity of the bone present. Note that we donot model diffusion in this equation as we assume that osteoblasts diffuse on a significantly lower levelthan the bio-active molecules. Of course, more than one cell type is present and responsible for bonegrowth. For simplicity we only include osteoblasts in this model, but an extension is easily feasible here. N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 5
Finally, the equation modelling bone growth (1.4) follows the same pattern as the one for osteoblastconcentration. In summary, our system of equations reads0 = div (cid:0) C ( ρ, σ, b ) ε ( u ) (cid:1) (mechanical equilibrium)(1.1) d t a i = div (cid:0) D i ( ρ ) ∇ a i (cid:1) + k ,i | ε ( u ) | δ c − k ,i a i (diffusion, generation, and decay of i =1 . . . N bio-molecules)(1.2) d t c = H ( a , . . . , a N , c, b ) (cid:18) − c − ρ (cid:19) (osteoblast generation)(1.3) d t b = K ( a , . . . , a N , c, b ) (cid:18) − b − ρ (cid:19) (bone regeneration driven by a, b and c ).(1.4)In the above system k , k ,i , k ,i ≥ i = 1 , . . . , N are constants that need to be determined fromexperiments, compare to the section 4 where we discuss certain choices. The functional relationships C , D i ( ρ ) , | · | δ , H and K are all required to satisfy certain technical assumptions that guarantee thewell-posedness of the above system. We discuss this in detail in section 2.Finally, we need to specify boundary conditions. For the elastic equilibrium equation we allow mixedboundary conditions including the limiting cases of a pure displacement boundary condition and a purestress boundary condition. As for the bio-active molecules we assume that these are in saturation, i.e., a ( t, x ) = 1 adjacent to bone and on the rest of the boundary of Ω we assume no-flux boundaryconditions. For the initial time-point we propose a i (0 , x ) = a i, = 0 inside of Ω. This choice reflects thescenario of a scaffold that is not preseeded with exogenous growth factors. However, different choices of a i, are admissible and allow the model to cover e.g., pre-seeding with osteoinductive factors. Finally, atthe initial time we assume that no osteoblasts and no regenerated bone are present inside the domainof computation. In formulas, it holds for all i = 1 , . . . , Na i (0 , x ) = 0 for all x ∈ Ω(1.5) a i ( t, x ) = 1 for all t ∈ I , x adjacent to bone(1.6) D ρi ∇ a i ( t, x ) · η = 0 for all t ∈ I , x not adjacent to bone(1.7) (cid:0) C ( ρ, σ, b ) ε ( u ( t, x )) (cid:1) · η = g N ( x ) on the Neumann boundary of Ω(1.8) u ( t, x ) = g D ( x ) on the Dirichlet boundary of Ω(1.9) c (0 , x ) = b (0 , x ) = 0 for all x ∈ Ω.(1.10)The model allows for a time dependent choice of the mechanical loading g D and g N . Due to the longregeneration time horizon of approximately 12 months, however, it is not expedient to resolve veryshort time-scales of, e.g., the mechanics of physical therapy. Instead, we consider suitably time-averagedloading conditions here.1.3. Concrete Examples.
We provide a number of possibilities for choosing the functional relation-ships C , D, H and K and boundary conditions for the mechanical equilibrium equation 1.1. For aneasy example of the elastic tensor that does not need to be derived by a complicated homogenizationprocedure we simply use the Voigt bound. If we denote by C b and C ρ the elastic tensors of maturedbone and intact PCL respectively (in their simplest form modelled as isotropic materials) we thuschoose C ( ρ, σ, b ) = b C b + ρσ C ρ . This is in accordance with [39] where the same idea was used in a model with only one spatial variable.Note that this C naturally is time-dependent as the quantities b and σ vary in time. While this examplemay serve as a first choice, one could also fix a concrete scaffold micro-structure, such as a gyroid design,and derive the explicit homogenized material properties (see, e.g., [3]). PATRICK DONDL, PATRINA S. P. POH, AND MARIUS ZEINHOFER
For the diffusivities D i ( ρ ) we propose a dependence on the scaffold density ρ , for example D i ( ρ ) = k i (1 − ρ ) Idwhere k i are constants that measure the diffusivity of the bio-active molecule a i without the presenceof the scaffold ρ . The term (1 − ρ ) accounts for reduced diffusivity for high PCL volume fractions. It isheuristically clear, yet interesting to note, that a too dense scaffold impairs bone regeneration. This isreflected in our model through the diffusivity above, since the amount of bioactive molecules is linkedto bone regeneration via the ODE (1.4). One could also imagine to derive the tensor D i ( ρ ) througha homogenization process which would then again reflect the choice of a specific micro-structure. Formathematical well-posedness reasons we are unable to allow the diffusivity D i ( ρ ) to depend on thebone density b . Furthermore, we also assume that D i ( ρ ) does not depend on time.Finally, we consider the functional relationships H and K inducing the production and proliferationof osteoblasts and bone. To be covered by our mathematical analysis, in the realization of H and K not more than two of the bio-active molecules should be multiplied. This is a technical mathematicalissue due to a possible lack of integrability. Compare also to Assumption 3.1 where we discuss thisissue rigorously. Consequently, we provide an example involving two bio-active molecules a and a .These can be assumed to have different production rates and half-lives. Then we set H ( a , a , c, b ) = H ( a , a , c ) = k a a (1 + k c )(1.11)hence bone growth only takes place when the full bio-environment, i.e., both molecules a and a arepresent. Furthermore the proliferation of osteoblasts is represented by the term (1 + k c ). Again k and k are some constants that need to be chosen in accordance with experiments.For K we propose a similar equation, modelling that bone growth takes place given the presence ofosteoblasts and a suitable biological environment, represented in the choice of K through the factor a . More precisely we set K ( a , a , c, b ) = K ( a , c ) = k a c. (1.12)Another choice for K reflecting that different bio-active molecules are responsible for different stagesof bone formation and maturation is possible. This makes the functional relationship dependent of b .We set K ( a, b ) = f ( b ) a c + f ( b ) a c. (1.13)Now, f can be chosen with support on small values of b , such that in this stage molecule a is drivingthe growth, and f with support on larger b , thus requiring a in later stages of regeneration. Weremark that empirically many different bio-molecules are observed and it is assumed that these arelinked to different biological processes [27].2. Mathematical Formulation
In this section we describe the mathematical setting in which we prove the existence of a solution tothe system of equations (1.1) – (1.4). We also state the assumptions the functional relationships C , D i , | · | δ H and K are required to satisfy.2.1. The Domain.
Fix a time interval I = [0 , T ] with T >
0. The spatial domain Ω ⊂ R n , with n = 1 , , ∂ Ω into a Dirichlet part and a Neumann part. For the elastic equation we write Γ eD and Γ eN for Dirichletand Neumann boundary respectively, here Γ eD = ∅ is allowed. For the diffusion equations we write Γ dD and Γ dN . To simplify notation we do not treat the case of different Dirichlet-Neumann partitions fordifferent diffusion equations, though this does not lead to further mathematical complications. Finallywe need to assume some regularity on Ω and the partition ∂ Ω = Γ dD ∪ Γ dN for the diffusion equations,namely the set Ω ∪ Γ dN needs to be Gr¨oger regular which is a concept introduced in [22], see also[24]. These regularity assumptions are tailored to provide a certain regularity of the solutions of the N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 7 diffusion equations which we discuss in detail in Appendix A. These assumptions are very general andcover the cases one wants to use in practice.2.2.
Admissible Data.
The admissible scaffold volume fractions ρ are given as(2.1) P := { ρ ∈ C (Ω) | c P ≤ ρ ( x ) ≤ C P } with some fixed constants 0 < c P < C P <
1, excluding unreasonable scaffold designs. To a scaffoldvolume fraction ρ ∈ P we assign the set W ρ of admissible cell and bone volume fractions, consistingof tuples of continuous functions in time and space(2.2) W ρ := { ( c, b ) ∈ C ( I × Ω) | ≤ c ( t, x ) , b ( t, x ) ≤ − ρ ( x ) } . The Elastic Equation.
We begin with the Hookean law C . It depends on the scaffold and bone,i.e., on ρ , σ and b and varies therefore in space and time. We assume that the map(2.3) W ρ → L ∞ ( I, L ∞ (Ω , L ( M s ))) with b (cid:55)→ ( t (cid:55)→ ( x (cid:55)→ C ( ρ, σ, b )( t, x )))is Lipschitz continuous with Lipschitz constant L C independent of ρ ∈ P . Remember that σ is a fixed ex-ponential decay. Here M s denotes the symmetric n × n matrices and L ( M s ) is the space of linear mapsfrom M s into itself, usually called the space of fourth order tensors. The space L ∞ ( I, L ∞ (Ω , L ( M s )))denotes a Bochner space, i.e., a Banach-space valued Lebesgue space, see, e.g., [18, 7]. In the follow-ing we will often omit the cumbersome notation of dependencies on x and t for C . Spelling out thedefinitions of the norms in (2.3) this Lipschitz continuity means that for all M ∈ M s it holds(2.4) | C ( ρ ( x ) , σ ( t ) , b ( t, x )) M − C ( ρ ( x ) , σ ( t ) , b ( t, x )) M | ≤ L C (cid:107) b − b (cid:107) C | M | for all ( c , b ), ( c , b ) ∈ W ρ and uniformly in ρ ∈ P and uniformly on the complement of a set ofmeasure zero in I × Ω. Furthermore we assume that there are constants 0 < c C < ∞ and 0 < C C < ∞ such that(2.5) sup ρ,c,b (cid:107) C ( ρ, σ, b ) (cid:107) L ∞ ( I,L ∞ (Ω , L ( M s ))) ≤ C C and inf ρ,c,b C ( ρ, σ, b ) M : M ≥ c C | M | where the supremum and infimum run over ρ ∈ P and b ∈ W ρ and A : B = tr AB T denotes thefull contraction of matrices. We now discuss the weak formulation of equation (1.1). Let ρ ∈ P and( c, b ) ∈ W ρ be some admissible functions. We first address the case where Γ eD has non-vanishingmeasure and comment on the pure Neumann problem later. The strong form − div (cid:16) C ( ρ, σ, b ) ε ( u ) (cid:17) = 0 in Ω , u | Γ eD = g eD , (cid:16) C ( ρ, σ, b ) ε ( u ) (cid:17) η | Γ eN = g eN encodes that at every point in time mechanical equilibrium is achieved, making the equation timedependent. The function space for the weak formulation is: L ( I, H , (Ω , R n )) with H , (Ω , R n ) beingthe Sobolev space of R n -valued, square integrable functions with square integrable derivatives, see forexample [8, 21, 1] for a detailed account of such spaces. If the context is clear, we will usually write H (Ω) instead of H , (Ω , R n ). The space of test functions is L ( I, H D e (Ω)), where H D e (Ω) is thesubspace of H (Ω) whose members vanish on Γ eD . For the Dirichlet boundary values we require g eD tobe in L ( I, H / (Γ eD , R n )), with H / (Γ), for some Γ ⊂ ∂ Ω, being the trace space of H (Ω), see forexample [1, 21]. The Neumann boundary values can be given as an element of L ( I, H / (Γ eN , R n ) (cid:48) ).Denoting by (cid:104)· , ·(cid:105) H / the dual pairing of H / (Γ eN , R n ) the weak formulation of (1.1) is (cid:90) I (cid:90) Ω C ( ρ, σ, b ) ε ( u ) : ε ( · ) dxdt = (cid:90) I (cid:104) g eN , ·(cid:105) H / dt in L ( I, H D e (Ω)) (cid:48) (2.6) u = g eD in L ( I, H / (Γ eD )) . The left hand side of (2.6) equation defines an operator T : L ( I, H (Ω)) → L ( I, H (Ω)) (cid:48) with T u = (cid:90) I (cid:90) Ω C ( ρ, σ, b ) ε ( u ) : ε ( · ) dxdt. PATRICK DONDL, PATRINA S. P. POH, AND MARIUS ZEINHOFER
Note that the isometry L ( I, H (Ω)) (cid:48) ˜= L ( I, H (Ω) (cid:48) ) implies that the equation (2.6) can be under-stood to hold almost everywhere in time, which is precisely what we want for our model. Furthermore,Korn’s inequality can be used to show that T is coercive, see [13]. The advantage of the abstract for-mulation is that it makes the Lax-Milgram Lemma applicable. Now we comment on the pure Neumannboundary value problem, i.e., the case Γ eN = ∂ Ω. We define the spaces W := ker( ε ) ⊂ H (Ω) and thequotient H (Ω) /W . Note that W consists of the functions of the form w ( x ) = Ax + b , where A is ananti-symmetric matrix and b ∈ R n , see for example [12]. For the pure Neumann problem consider theoperator T : L ( I, H (Ω) /W ) → L ( I, H (Ω) /W ) (cid:48) using the induced map ˆ ε : H (Ω) /W → L (Ω , M s ) in its definition T ( u ) = (cid:90) I (cid:90) Ω C ( ρ, σ, b )ˆ ε ( u ) : ˆ ε ( · ) dxdt. The codomain of this operator is L ( I, H (Ω) /W ) (cid:48) ˜= L ( I, ( H (Ω) /W ) (cid:48) ), which encodes a compat-ibility condition. We assume that our Neumann boundary condition is given as a function g eN ∈ L ( I, H / ( ∂ Ω) (cid:48) ) that satisfies almost everywhere in I (2.7) (cid:104) g eN ( t ) , ·(cid:105) H / = 0 for all w ∈ W. This guarantees that (cid:90) I (cid:104) g eN , · (cid:105) H / dt ∈ L ( I, H (Ω) /W ) (cid:48) is an admissible right hand side. The pure Neumann problem consists then of finding u ∈ L ( I, H (Ω) /W )such that (cid:90) I (cid:90) Ω C ( ρ, σ, b )ˆ ε ( u ) : ˆ ε ( · ) dxdt = (cid:90) I (cid:104) g eN , · (cid:105) H / dt ∈ L ( I, H (Ω) /W ) (cid:48) , Finally, let us remark that one can treat the Dirichlet, the Neumann and the mixed boundary valueproblem at once by always passing to the quotient H D e (Ω) /W . In the case of a proper Dirichletboundary condition we then have W ∩ H D e (Ω) = { } , which implies H D e (Ω) /W = H D e (Ω), hencerecovers the Dirichlet or mixed case, and if Γ eN = ∂ Ω we retrieve the pure Neumann case.2.4.
Diffusion Equations.
Before we state the weak formulation of the diffusion equations, for thereader’s convenience, we recall the concept of the time derivative we are using – namely a regularBanach space valued distribution with a dense embedding j ∈ L ( X, X (cid:48) ) just as in [7]. Let ( i, X, H )be a Gelfand triple, i.e., X is a Banach space, H is a Hilbert space and i ∈ L ( X, H ) has dense range.Then we set j to be j = i (cid:48) ◦ R ◦ i where R : H → H (cid:48) is the Riesz isometry and i (cid:48) denotes the Banachspace adjoint of i . We say a function a ∈ L ( I, X ) possesses a time derivative d t a ∈ L ( I, X (cid:48) ) if itholds (cid:90) I ( j ◦ a )( t ) ∂ t ϕ ( t ) dt = − (cid:90) I d t a ( t ) ϕ ( t ) dt ∀ ϕ ∈ D ( I ) . The integrals are X (cid:48) valued Bochner integrals and we set D ( I ) := C ∞ c ( I ) as usual. This is used todefine a generalized Sobolev space built on the triple ( i, X, H ) as H , , ( I, X, X (cid:48) ) = { a ∈ L ( I, X ) | d t a ∈ L ( I, X (cid:48) ) } . See in [7, Chapter II, section 5] for more information. We only remark that functions in this Sobolevspace have representatives in C ( I, H ), hence initial value problems can be formulated.To get to our concrete diffusion equations we let ρ ∈ P , ( c, b ) ∈ W ρ and, depending on the boundaryconditions for the elastic equation, u ∈ L ( I, H (Ω)) or u ∈ L ( I, H (Ω) /W ) be some fixed functions.In order to work with homogeneous Dirichlet boundary conditions in space we write a i ( t ) = ˜ a i ( t ) + 1 with ˜ a i ( t ) ∈ H D d (Ω) for i = 1 , . . . , N. N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 9
Here H D d (Ω) denotes the subspace of H (Ω) with vanishing trace on Γ dD . We can thus seek ˜ a i in thespace H , , ( I, H D d (Ω) , H D d (Ω) (cid:48) ) built around the triple (id | H Dd , H D d , L ) satisfying the equation (cid:90) (cid:104) d t ˜ a i , ·(cid:105) H Dd + (cid:90) (cid:90) D ρi ∇ ˜ a i ∇ · + k i ˜ a i · dxdt = (cid:90) (cid:90) ( k i | ε ( u ) | δ c − k i ) · dxdt (2.8) ˜ a i (0) = − . (2.9)The first equation is an equality in the space L ( I, H D d (Ω)) (cid:48) , i.e., it is required to hold when testedwith all members of L ( I, H D d (Ω)). In the second equation, the initial conditions is an equality in thespace L (Ω). For every i = 1 , . . . , N we have different constants k i and k i and also different diffusivities D ρi . Note that the quantity | ε ( u ) | d elta is well defined, even though the solution of the elastic equationis only unique up to rigid body motions. We assume furthermore that the D ρi are time-independent,measurable, essentially bounded and coercive, precisely D ρi ∈ L ∞ (Ω , M s )(2.10) (cid:104) D ρi ξ, ξ (cid:105) ≥ c D | ξ | ∀ ξ ∈ R n (2.11)where M s again denotes the symmetric n × n matrices and the inequality in (2.11) is to be understooduniformly in x ∈ Ω, ρ ∈ P and i = 1 , . . . , N . Finally the function | · | δ : R n × n → [0 , ∞ ) is required toto be globally Lipschitz and to satisfy an estimate of the form(2.12) | A | δ ≤ C | A | + C for all A ∈ R n × n where C , C > | A | denotes the Euclidean norm of a matrix.2.5. Ordinary Differential Equations.
We treat the ordinary differential equations in the vectorvalued sense and focus here on the cell equation (1.3), the bone equation (1.4) being treated analogously.For each x ∈ Ω, we thus seek a function c x satisfying the ODE c (cid:48) x ( t ) = H ( a ( t, x ) , . . . , a N ( t, x ) , c x ( t ) , b ( t, x )) (cid:18) c x ( t )1 − ρ ( x ) (cid:19) with c x (0) = 0. If there is a solution for all x ∈ Ω we obtain a function c in time and space, i.e., c : I × Ω → R with c ( t, x ) := c x ( t ). As H ( a , . . . , a N , c, b ) can not generally assumed to be continuous,a reasonable space to work in is W ,p ( I, X ) = { c ∈ L p ( I, X ) | d t c ∈ L p ( I, X ) } , similar to the space for the diffusion equation, but without the identification j : X (cid:44) → X (cid:48) . An existenceand uniqueness result in this setting can be found in the appendix, see Theorem B.2.In our concrete case we choose X = C (Ω), p = 2, so for fixed ρ ∈ P and a = ( a , . . . , a N ) ∈ H , , ( I, H (Ω) , H D d (Ω) (cid:48) ) N we seek c ∈ W , ( I, C (Ω)) satisfying(2.13) d t c = H ( a , . . . , a N , c, b ) (cid:18) − c − ρ (cid:19) with c (0) = 0 . We assume that H is a Nemytskii operator induced by a function which we again denote by H ,(2.14) H : R N +2 → R with ( a , . . . , a N ) = a (cid:55)→ H ( a )such that H ( a, c, b ) ≥ a , . . . , a N , b, c ≥
0. Furthermore we assume that H is locallyLipschitz continuous. Note that by some abuse of notation we denote by a , b and c both a function ina Sobolev space and a vector in Euclidean space.For the bone ODE we work in the same space and seek b ∈ W ,q ( I, C (Ω)) satisfying(2.15) d t b = K ( a , . . . , a N , c, b ) (cid:18) − b − ρ (cid:19) with b (0) = 0 . We assume the functional relationship K is induced by K : R N +2 → R with ( a, b, c ) = ( a , . . . , a N , b, c ) (cid:55)→ K ( a, b, c ) that satisfies K ( a , . . . , a N , b, c ) ≥ a , . . . , a N , b, c ≥ K is locally Lipschitz continuousas a map K : R N +2 → R . Finally, we need another assumption on H and K that is connected to theintegrability and the regularity properties of the solutions to the diffusion equations, see assumption3.1. We summarize our setting. Assumption 2.1.
We assume domain regularity as discussed in subsection 2.1, define the admissiblescaffold densities P in (2.1) and the set W ρ in (2.2). The material tensor C satisfies (2.3) and (2.5) andadmissible boundary conditions for the elastic equation are given in (2.6) and (2.7). For the diffusion weassume (2.10) and (2.11) and | · | δ must satisfy (2.12). The functional relationships H and K need to belocally Lipschitz, preserve positivity and satisfy the technical assumption 3.1 concerning integrability.3. Existence and Uniqueness
In this section we will prove that there exists a unique solution to the system (1.1)–(1.4) in the weaksense, i.e., there are functions u ∗ = ˜ u ∗ + u g eD with ˜ u ∗ ∈ L ( I, H D e (Ω) /W ) and u g eD | Γ eD = g eD , a ∗ = ˜ a ∗ +1with ˜ a ∗ ∈ H ( I, H D d (Ω) , H D d (Ω) (cid:48) ), c ∗ ∈ W ,p ( I, C (Ω)) and b ∗ ∈ W ,q ( I, C (Ω)) satisfying (cid:90) I (cid:90) Ω C ( ρ, σ, b ∗ )ˆ ε (˜ u ∗ + u g eD ) : ˆ ε ( · ) dxdt = (cid:90) I (cid:104) g eN , · (cid:105) H / (Γ eN ) dt (3.1) (cid:90) (cid:104) d t ˜ a ∗ i , ·(cid:105) + (cid:90) (cid:90) D ρi ∇ ˜ a ∗ i ∇ · + k i ˜ a ∗ i · dxdt = (cid:90) (cid:90) ( k i | ε ( u ∗ ) | δ c ∗ − k i ) · dxdt (3.2) ˜ a ∗ i (0) = − , with i = 1 , . . . , N, (3.3) d t c ∗ = H ( a ∗ , . . . , a ∗ N , c ∗ , b ∗ ) (cid:18) − c ∗ − ρ (cid:19) with c ∗ (0) = 0 , (3.4) d t b ∗ = K ( a ∗ , . . . , a ∗ N , c ∗ , b ∗ ) (cid:18) − b ∗ − ρ (cid:19) with b ∗ (0) = 0 . (3.5)The proof of this result relies essentially on the elementary fixed point theorem of Banach which wewill employ for the complete metric space W ρ . The strategy is to fix ρ ∈ P , then start with somearbitrary admissible functions ( c, b ) ∈ W ρ and to solve the equations successively. More precisely, theelastic equation will yield u = u ( c, b ), the diffusion equations a i = a i ( c, u ), the cell equation will besolved with data a i and b yielding an updated cell function c = c ( a i , b ) and finally the bone equationwill be solved with data a i and c to get an updated bone function b = b ( a i , c ). This procedure givesrise to an operator I which we will refer to as the iteration operator, formally I : W ρ → W ρ with ( c, b ) (cid:55)→ ( c, b ) . It is easy to see that all possible solutions to (3.1)–(3.5) correspond to all possible fixed-points of I . Thecrucial part of the proof consists of establishing regularity for the solutions of the diffusion equations,see also Appendix A for a discussion of results known in the literature serving our purpose.Finally, the whole strategy discussed above does only work on a short time interval I = [0 , T ], i.e., T small enough. However, by a continuation argument we can afterwards extend this solution to spanany finite time interval. We will need a technical assumption on the ODEs in connection with theiteration operator I . This is due to the fact that we cannot guarantee an L ∞ ( I × Ω) bound on thesolutions to the diffusion equations. See also Remark 3.3 on when the following assumption holds.
Assumption 3.1.
Let ρ ∈ P and ( c, b ) ∈ W ρ and denote by u ∈ L ( I, H (Ω) /W ), a ∈ L ( I, C (Ω) N )and c ∈ C ( I × Ω) the functions produced by solving the equations successively as in the definitionof the iteration operator I . The existence and regularity of these solutions is discussed in the maintheorem. Assume there exists p ∈ [1 , ∞ ] such that for every bounded set B ⊂ C (Ω) there are functions m HB ∈ L p ( I ) and L HB ∈ L ( I ) such that it holds (cid:107) H ( a ( t ) , ˜ c, b ( t ) (cid:107) C (Ω) ≤ m HB ( t ) for all ˜ c ∈ B, a.e. in I, (3.6) (cid:13)(cid:13) H ( a ( t ) , ˜ c, b ( t )) − H ( a ( t ) , ˜˜ c, b ( t )) (cid:13)(cid:13) C (Ω) ≤ L HB ( t ) for all ˜ c, ˜˜ c ∈ B, a.e. in I. (3.7) N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 11
Furthermore assume that there are functions λ H and β H in L ( I ) and α H ∈ L ( I ) such that we canestimate, independently of the choice of ( c, b ) ∈ W ρ (and consequently a and c ), (cid:107) H ( a ( t ) , c ( t ) , b ( t )) (cid:107) C (Ω) ≤ λ H ( t ) . (3.8)Additionally, uniformly for any ( c , b ) , ( c , b ) ∈ W ρ and corresponding a , a , c , c , we have (cid:13)(cid:13) H ( a ( s ) , b ( s ) , c ( s )) − H ( a ( s ) , b ( s ) , c ( s )) (cid:13)(cid:13) C (Ω) ≤ β H ( t ) (cid:107) c ( s ) − c ( s ) (cid:107) C (Ω) + α H ( t ) (cid:104) (cid:13)(cid:13) a ( s ) − a ( s ) (cid:13)(cid:13) C (Ω) N + (cid:107) b ( s ) − b ( s ) (cid:107) C (Ω) (cid:105) . (3.9)For K we assume analogous properties, i.e., there is q ∈ [1 , ∞ ] such that for every bounded B ⊂ C (Ω)there are m KB ∈ L q ( I ) and L KB ∈ L ( I ) as well as λ K , β K ∈ L ( I ) and α K ∈ L ( I ) satisfying estimatesas above. Theorem 3.2 (Existence & Uniqueness) . Let ρ ∈ P be fixed and let the Assumptions 2.1 and 3.1 befulfilled. Then there exist unique functions u ∗ = ˜ u ∗ + u g eD with ˜ u ∗ ∈ L ( I, H D e (Ω) /W ) and u g eD | Γ eD = g eD , a ∗ = ˜ a ∗ + 1 with ˜ a ∗ ∈ H ( I, H D d (Ω) , H D d (Ω) (cid:48) ) , c ∗ ∈ W ,p ( I, C (Ω)) and b ∗ ∈ W ,q ( I, C (Ω)) solving the system (3.1) – (3.5) .Proof. We need to establish the contraction and self-mapping property of I . Let us thus fix two tuples( c , b ) and ( c , b ) ∈ W ρ . We aim to show an estimate of the form (cid:107)I ( c , b ) − I ( c , b ) (cid:107) C ( I × Ω) = (cid:107) c − c (cid:107) C ( I × Ω) + (cid:13)(cid:13) b − b (cid:13)(cid:13) C ( I × Ω) ≤ C ( I ) (cid:0) (cid:107) c − c (cid:107) C ( I × Ω) + (cid:107) b − b (cid:107) C ( I × Ω) (cid:1) , where C ( I ) → | I | →
0, making I the desired self-mapping for T small enough. The Elastic Equation.
We will treat a pure Neumann and a mixed boundary value problem simul-taneously. We endow the space H D e (Ω) /W with the norm (cid:107) ε ( · ) (cid:107) L (Ω) , which by Korn’s inequality isequivalent to the natural one on H D e (Ω) /W , see for example [13]. By definition of g eD , there is a func-tion u g eD ∈ L ( I, H (Ω) /W ) such that u g eD | Γ eD = g eD . In the weak formulation of the elastic equationwe seek ˜ u i ∈ L ( I, H D e (Ω) /W ) satisfying(3.10) (cid:90) I (cid:90) Ω C ( ρ, σ, b i )ˆ ε (˜ u i + u g eD ) : ˆ ε ( · ) dxdt = (cid:90) I (cid:104) g eN , · (cid:105) H / (Γ eN ) dt in the space L ( I, H D e (Ω) /W ) (cid:48) . Then u i := ˜ u i + u g eD is the solution we are interested in. Note that ifΓ eD has vanishing measure, we can choose u g eD = 0 and H D e (Ω) /W = H (Ω) /W . On the other hand,if Γ eD has positive measure, then W ∩ H D e (Ω) = 0 and H D e (Ω) /W = H D e (Ω). The equation (3.10)leads to the operators T b i : L ( I, H D e (Ω) /W ) → L ( I, H D e (Ω) /W ) (cid:48) with T b i v = (cid:90) I (cid:90) Ω C ( ρ, σ, b i )ˆ ε ( v ) : ˆ ε ( · ) dxdt and right hand sides f b i = (cid:90) I (cid:104) g eN , · (cid:105) H / (Γ eN ) dt (cid:124) (cid:123)(cid:122) (cid:125) =: f N − (cid:90) I (cid:90) Ω C ( ρ, σ, b i )ˆ ε ( u g eD ) : ˆ ε ( · ) dxdt (cid:124) (cid:123)(cid:122) (cid:125) =: f Dbi . By our assumption (2.5) and Korn’s inequality the operators T b i are coercive with coercivity con-stant c C . Applying the Lax-Milgram Lemma we find that there are unique solutions ˜ u and ˜ u ∈ L ( I, H D e (Ω) /W ) to T b i ˜ u i = f b i . By the duality L ( I, H (Ω) /W ) (cid:48) = L ( I, ( H (Ω) /W ) (cid:48) ) we know that almost everywhere in I the function u i ( t ) = ˜ u i ( t ) + u g eD ( t ) satisfies (cid:90) Ω C ( ρ, σ, b i )( t )ˆ ε (˜ u i )( t ) : ˆ ε ( · ) dx = (cid:104) g eN ( t ) (cid:105) H / (Γ eN ) − (cid:90) Ω C ( ρ, σ, b i )ˆ ε ( u g eD ( t )) : ˆ ε ( · ) dx in the space H (Ω) /W . Using Lax-Milgram again we get using the boundedness and coercivity con-stants from (2.5) (cid:107) u i ( t ) (cid:107) H (Ω) /W ≤ c − C (cid:104) (cid:107) g eN ( t ) (cid:107) H / (Γ eN ) (cid:48) + C C (cid:13)(cid:13) u g eD ( t ) (cid:13)(cid:13) H (Ω) /W (cid:105) . As the above estimate is independent of ρ, c i and b i it holdssup ρ,c,b (cid:107) u ( ρ, b ) (cid:107) L ( I,H (Ω) /W ) ≤ C ( I ) , (3.11)where u ( ρ, c, b ) denotes the solution of the elastic problem to the data ρ ∈ P and ( c, b ) ∈ W ρ . To showthat C ( I ) tends to zero with | I | → u − u . We claim that(3.12) (cid:107) u − u (cid:107) L ( I,H (Ω) /W ) ≤ C ( I ) (cid:107) b − b (cid:107) C ( I × Ω) where again C ( I ) → | I | →
0. To establish this, note that ˜ u − ˜ u = u − u and compute f b − f b = f Db − f Db = T b ˜ u − T b ˜ u = T b (˜ u − ˜ u ) + T b ˜ u − T b ˜ u − T b ˜ u . Hence T b ( u − u ) = ( T b ˜ u − T b ˜ u ) + ( f Db − f Db ) and using (cid:13)(cid:13) T − b (cid:13)(cid:13) ≤ c − C we find (cid:107) u − u (cid:107) L ( I,H (Ω) /W ) ≤ c − C (cid:107)T b ˜ u − T b ˜ u (cid:107) L ( I,H (Ω) /W ) (cid:48) + c − C (cid:13)(cid:13) f Db − f Db (cid:13)(cid:13) L ( I,H (Ω) /W ) (cid:48) . We estimate the terms of the right hand side using the Lipschitz continuity of C which we assumed in(2.4), combining it with (3.11) to find (cid:107)T b ˜ u − T b ˜ u (cid:107) L ( I,H (Ω) /W ) (cid:48) ≤ L C (cid:107) b − b (cid:107) C ( I × Ω) (cid:107) ˜ u (cid:107) L ( I,H (Ω) /W ) ≤ L C C ( I ) (cid:107) b − b (cid:107) C ( I × Ω) = C ( I ) (cid:107) b − b (cid:107) C ( I × Ω) and (cid:13)(cid:13) f Db − f Db (cid:13)(cid:13) L ( I,H (Ω) /W ) (cid:48) ≤ L C (cid:13)(cid:13) u g eD (cid:13)(cid:13) L ( I,H (Ω) /W ) (cid:107) b − b (cid:107) C ( I × Ω) = C ( I ) (cid:107) b − b (cid:107) C ( I × Ω) . The Diffusion Equations.
Given the functions c i , u i with i = 1 , ρ , we turn to the diffusionequations. We seek functions a i = ˜ a i + 1 where the ˜ a i are members of H ( I, H D d (Ω) , H D d (Ω) (cid:48) ) N , thatmeans a i = ( a i , . . . , a iN ) = (˜ a i + 1 , . . . , ˜ a iN + 1), i = 1 ,
2, denoting the components of a i with lowerindices. For j = 1 , . . . , N the ˜ a ij are sought to satisfy the following equation in L ( I, H D d (Ω)) (cid:48) (cid:90) I (cid:104) d t ˜ a ij , · (cid:105) H D dt (cid:124) (cid:123)(cid:122) (cid:125) =: d t ˜ a ij + (cid:90) (cid:90) D ρj ∇ ˜ a ij ∇ · + k j ˜ a ij · dxdt (cid:124) (cid:123)(cid:122) (cid:125) =: M j ( ρ )˜ a ij = (cid:90) (cid:90) ( k j | ε ( u i ) | c i − k j ) · dxdt (cid:124) (cid:123)(cid:122) (cid:125) =: f jui,ci and initial value ˜ a j (0) = − L (Ω). The operators( d t + M j ( ρ ) , ev ) : H ( I, H D d (Ω) , H D d (Ω) (cid:48) ) → L ( I, H D d (Ω)) (cid:48) × L (Ω)are linear homeomorphisms, see for example [20] for a proof, which essentially relies on the coercivity of M j ( ρ ). This explains why we assumed (2.11) and hence we can guarantee the existence of the ˜ a ij . We N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 13 state two important properties of the solutions a ij and their differences a j − a j , to which references orproofs can be found in the Appendix A. The first is a lower pointwise bound, it holds for j = 1 , . . . , N and i = 1 , ≤ a ij ( t, x ) = a ij ( t, x ) almost everywhere in I × Ω . (3.13)This is due to the positivity of the right hand sides f ju i ,c i . Secondly, we look at the equations satisfiedby the differences a j − a j . These equations possess right hand sides f ju ,c − f ju ,c in L ( I, L (Ω))and with ( a j − a j )(0) = 0 smooth initial conditions. Then, using regularity for mixed boundary valueproblems, see Theorem A.1, there is α > (cid:13)(cid:13) a j − a j (cid:13)(cid:13) L ( I,C α (Ω)) ≤ C (cid:13)(cid:13) f ju ,c − f ju ,c (cid:13)(cid:13) L ( I,L (Ω)) . (3.14)The constant C is uniform in the data ρ ∈ P , ( c, b ) ∈ W ρ and u ( ρ ). We claim now that we get thefollowing estimate for the difference a − a (cid:13)(cid:13) a − a (cid:13)(cid:13) L ( I,C α (Ω)) N ≤ C (cid:16) (cid:107) c − c (cid:107) C ( I × Ω) + (cid:107) u − u (cid:107) L ( I,H /W ) (cid:17) (3.15)with C not blowing up as | I | →
0. This estimate is obtained, using (3.14) and estimating the difference f ju ,c − f ju ,c . It holds f jc ,u − f jc ,u = k j | ε ( u ) | δ ( c − c ) + k j ( | ε ( u ) | δ − | ε ( u ) | δ ) c . Using the fact that c takes values in the unit interval and the assumptions on | · | δ , see 2.12, it follows (cid:13)(cid:13) f jc ,u − f jc ,u (cid:13)(cid:13) L ( I,L (Ω)) ≤ C (cid:0) (cid:107) ε ( u ) (cid:107) L ( I,L (Ω)) + 1 (cid:1) (cid:107) c − c (cid:107) C ( I × Ω) + C (cid:107) ε ( u − u ) (cid:107) L ( I,L (Ω)) . Invoking (3.11) we know that (cid:107) ε ( u ) (cid:107) L ( I,L (Ω)) is bounded uniformly in the data ρ ∈ P and ( c, b ) ∈ W ρ .Combining this with the identity (cid:107) ε ( u − u ) (cid:107) L ( I,L (Ω)) = (cid:107) u − u (cid:107) L ( I,H /W ) we conclude. The Cell ODE.
We turn now to the Cell ODE and solve this equation twice, once with data ρ, a , . . . , a N and b , producing a function c , and once with ρ, a , . . . , a N and b yielding c . The solutions c and c are members of the space W ,p ( I, C (Ω)) and consequently of C ( I × Ω) satisfying 0 ≤ c ( t, x ) ≤ − ρ ( x )solving the ODE d t c i = H ( a i , . . . , a iN , b i , c i ) (cid:18) − c i − ρ (cid:19) with c i (0) = 0 . These facts are proven as Lemma B.6 in the Appendix. Our goal is to estimate the difference c − c and we claim that it holds(3.16) (cid:107) c − c (cid:107) C ( I × Ω) ≤ C ( I ) (cid:16) (cid:13)(cid:13) a − a (cid:13)(cid:13) L ( I,C α (Ω) N ) + (cid:107) b − b (cid:107) C ( I × Ω) (cid:17) where C ( I ) tends to zero with | I | →
0. To prove the estimate 3.16 we use the fundamental theorem ofthe space W ,p ( I, C (Ω)) and write c i ( t ) = (cid:90) t H ( a i ( s ) , b i ( s ) , c i ( s )) (cid:18) − c i ( s )1 − ρ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) =: γ i ( s ) ds. By Assumption 3.1 the expression γ ( s ) − γ ( s ) can be estimated (cid:107) γ ( s ) − γ ( s ) (cid:107) C (Ω) ≤ (cid:13)(cid:13) H ( a , b , c )( s ) − H ( a , b , c )( s ) (cid:13)(cid:13) C (Ω) (cid:13)(cid:13)(cid:13)(cid:13) − c ( s )1 − ρ (cid:13)(cid:13)(cid:13)(cid:13) C (Ω) + (cid:13)(cid:13) H ( a , b , c )( s ) (cid:13)(cid:13) C (Ω) (cid:13)(cid:13) (1 − ρ ) − (cid:13)(cid:13) (cid:107) c ( s ) − c ( s ) (cid:107) C (Ω) ≤ Cα H ( s ) (cid:104) (cid:13)(cid:13) a ( s ) − a ( s ) (cid:13)(cid:13) C (Ω) N + (cid:107) b ( s ) − b ( s ) (cid:107) C (Ω) (cid:105) + C ( β H ( s ) + λ H ( s )) (cid:107) c ( s ) − c ( s ) (cid:107) C (Ω) . Now we can apply Gr¨onwall’s and H¨older’s inequality, using that α H ∈ L ( I ) and β H , λ H ∈ L ( I ) toobtain (cid:107) c ( t ) − c ( t ) (cid:107) C (Ω) ≤ C (cid:90) I α H ( s ) (cid:104) (cid:13)(cid:13) a ( s ) − a ( s ) (cid:13)(cid:13) C (Ω) N + (cid:107) b ( s ) − b ( s ) (cid:107) C (Ω) (cid:105) ≤ C (cid:13)(cid:13) α H (cid:13)(cid:13) L ( I ) (cid:104) (cid:107) a − a (cid:107) L ( I,C (Ω) N ) + (cid:107) b − b (cid:107) L ( I,C (Ω)) (cid:105) ≤ C ( I ) (cid:104) (cid:13)(cid:13) a − a (cid:13)(cid:13) L ( I,C α (Ω)) + (cid:107) b − b (cid:107) C ( I × Ω) (cid:105) . Here we used that (cid:13)(cid:13) α H (cid:13)(cid:13) L ( I ) → | I | →
0, which follows from Lebesgue’s dominated convergencetheorem. As the right side of the estimate is independent of t ∈ I , this shows that (3.16) holds. The Bone ODE.
Finally we treat the Bone ODE. Again we solve it twice, with data a i , . . . , a iN , c i and ρ producing b i with i = 1 ,
2. The functions b & b are members of W ,q ( I, C (Ω)) and conse-quently of C ( I × Ω) satisfying 0 ≤ b i ( t, x ) ≤ − ρ ( x ) and d t b i = K ( a i , . . . a iN , b i , c i ) (cid:18) − b i − ρ (cid:19) This means that b i ∈ W ρ , hence making the iteration map I a self mapping. All these properties areestablished as in the case of the Cell ODE. Repeating our computations for c − c we find (cid:13)(cid:13) b − b (cid:13)(cid:13) C ( I × Ω) ≤ C ( I ) (cid:16) (cid:13)(cid:13) a − a (cid:13)(cid:13) L ( I,C α (Ω) N ) + (cid:107) c − c (cid:107) C ( I × Ω) (cid:17) . (3.17)and the constant C ( I ) tends to zero as | I | → Contraction Property of I . We collect all estimates to see that I is a contractive self-mapping for | I | small enough. Use (3.17), (3.16), (3.15), and (3.12) to conclude (cid:13)(cid:13) b − b (cid:13)(cid:13) C ( I × Ω) ≤ C ( I ) (cid:16) (cid:13)(cid:13) a − a (cid:13)(cid:13) L ( I,C α (Ω) N ) + (cid:107) c − c (cid:107) C ( I × Ω) (cid:17) ≤ C ( I ) (cid:16) (cid:13)(cid:13) u − u (cid:13)(cid:13) L ( I,H (Ω) /W ) + (cid:107) c − c (cid:107) C ( I × Ω) (cid:17) ≤ C ( I ) (cid:16) (cid:107) b − b (cid:107) C ( I × Ω) + (cid:107) c − c (cid:107) C ( I × Ω) (cid:17) . and the estimate for (cid:107) c − c (cid:107) C ( I × Ω) works identically. Consequently, it holds (cid:107)I ( c , b ) − I ( c , b ) (cid:107) C ( I × Ω) ≤ C ( I ) (cid:104) (cid:107) c − c (cid:107) C ( I × Ω) + (cid:107) b − b (cid:107) C ( I × Ω) (cid:105) As C ( I ) → | I | →
0, the contraction map principle implies that I : ( c, b ) (cid:55)→ ( c, b ) possesses aunique fix point for | I | small enough. Long-Time Existence.
We established the existence of a solution ( u ∗ , a ∗ , c ∗ , b ∗ ) on an interval [0 , T ]where T > I a contraction. Now we use the well defined functions c ∗ ( T, · )and b ∗ ( T, · ) ∈ C (Ω) as initial data for the ODEs and as a ∗ ∈ C ([0 , T ] , L (Ω) N ) the function a ∗ ( T ) ∈ L (Ω) N serves as start value for the diffusion equations. Repeating the computations we findthat there exists a unique solution ( u ∗∗ , a ∗∗ , c ∗∗ , b ∗∗ ) to the system on the interval [ T − ε, T − ε ] for N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 15 some small ε >
0. On the overlap [ T − ε, T ] the solutions ( u ∗∗ , a ∗∗ , c ∗∗ , b ∗∗ ) and ( u ∗ , a ∗ , c ∗ , b ∗ ) agreeand thus we found the unique solution on the interval [0 , T − ε ]. As T does not depend on the initialvalues of neither a ∗ , c ∗ nor b ∗ this iterates to span every finite time interval. (cid:3) Remark . We discuss the validity of Assumption 3.1, which is given in an implicit form. We treatroughly two cases. Either, only treating pure Dirichlet problems for the diffusion equations or assumingstrong properties for | · | δ , one can establish an L ∞ ( I × Ω) bound on the solutions of the diffusionequations and can then apply Lemma 3.4, or one is allowed to only multiply at most two differentcomponents of a in order not to violate the integrability. Furthermore we will always assume that weare in the setting of section 2.(i) Assume that | · | δ : R n × n → R is a bounded function. Then the solution to the diffusionequations lie in L ∞ ( I × Ω) with a bound on the uniform norm not depending on ρ ∈ P and b ∈ W ρ . In view of Lemma 3.4 this implies that Assumption 3.1 holds.(ii) Assume that we consider a pure Dirichlet problem for the diffusion equations and that theDirichlet data on the parabolic boundary lies in the space L ∞ ( I, L ∞ ( ∂ Ω)). Theorem 7.1 andCorollary 7.1 in [29] show that the solutions of the diffusion equations are members of L ∞ ( I × Ω)with a uniform bound on their norms. Here one crucially needs the Dirichlet information onthe parabolic boundary, thus the assumptions. We currently do not know whether a similarresult is available in the case of mixed boundary conditions.(iii) If we do not assume anything besides the setting of section 2, in the choice of H and K notmore than two of the components of a should be multiplied. In particular the choice | · | δ = | · | is covered. A concrete example is given in section 4. Lemma 3.4.
Let Assumption 2.1 hold and assume that for any choice of ρ and b ∈ W ρ the function a ∈ H ( I, H D d (Ω) , H D d (Ω) (cid:48) ) N produced by the iteration operator I is a member of L ∞ ( I × Ω) N with abound on the L ∞ ( I × Ω) norm that is uniform in ρ ∈ P and b ∈ W ρ . Then Assumption 3.1 is satisfied.Proof. For a fixed bounded set B ⊂ C (Ω) the following subset of R N +2 { ( a ( t, x ) , ˜ c ( t, x ) , b ( t, x )) | ˜ c ∈ B, ( t, x ) ∈ I × Ω } is relatively compact. By the continuity of H we can choose m HB to be a constant function, i.e., amember of L ∞ ( I ). Using the Lipschitz continuity of H on the set defined above, we are able toestablish property (3.7). Now let b and b ∈ W ρ be given and correspondingly a , a , c and c . Usingthat the set { ( a i ( t, x ) , c i ( t, x ) , b i ( t, x )) | i = 1 , t, x ) ∈ I × Ω } is bounded independently of ρ, b , b etc. we may use again the above reasoning and obtain that (3.1)holds with a constant function λ H and (3.1) with constant functions β H and α H . The remainingrequirements in assumption 3.1 are satisfied likewise. (cid:3) Numerical Experiments
In [15] porous PCL scaffolds with a periodic honeycomb structure and 87% porosity were used as atreatment strategy for 30mm tibial defects in an ovine model. This experiment was conducted in twogroups, one preseeding the scaffold with a special bio-active molecule (BMP) and the second groupwithout such preseeding. Here, we aim to numerically recreate the experiment without preseeding,using a concrete instance of our computational model.As usual, the experimental setup in [15] includes the use of a so-called fixateur – a titanium or steelplate that is fixed to the bone surrounding the defect site using screws. This fixateur is used to provideadditional mechanical stability. We include this device in a simplistic manner in our simulations,neglecting the effect of screws. From a modeling perspective, the fixateur acts as a stress shield on oneside of the defect and thus influences bone growth significantly.
As a concrete instance of our model we use two bioactive molecules and consider the following systemof equations 0 = div (cid:16) C ( ρ, σ, b ) ε ( u ) (cid:17) d t a = div (cid:16) D ( ρ ) ∇ a (cid:17) + k , | ε ( u ) | c − k , a d t a = div (cid:16) D ( ρ ) ∇ a (cid:17) + k , | ε ( u ) | c − k , a d t c = k a a (1 + k c ) (cid:18) − c − ρ (cid:19) d t b = k a c (cid:18) − b − ρ (cid:19) . We use mixed boundary values for the elastic equilibrium equation, with a surface traction stemmingfrom a force of 0 . . a ( t, x ) = a ( t, x ) = 1, adjacent to boneand otherwise we assume a non-flux boundary condition. Osteoblast and bone density is set to zero atthe initial time-point. Note that the concrete choice of boundary conditions here should be considereda proof of concept. Further, more detailed numerical studies are forthcoming.4.1. Model Parameters.
We report the choices for the constants and functional relationships intable 1. Some comments are in order.(a) In a healthy individual, given appropriate clinical interventions, bone defects should be com-pletely bridged with low to medium weight-bearing capacity after 6 months, see [52]. Thebone remodelling process to follow can take 3 to 5 years until the full function of the bone isrestored. We therefore consider a time span of 12 months for our model, which we identifiedas the critical phase for scaffold mediated bone healing.(b) The PCL decay parameter, k , is based on the experimental studies in [37], which shows thatafter one year 30% of the molecular mass remains.(c) The surface traction is set to 0 . . . We propose to view this time-constant surface traction as an averaged maximal stress.Furthermore we assume that due to the injury this averaged maximal stress is considerablylower than what is to be expected in a healthy individual, where we set it to 0 . . k ,i , k ,i , i = 1 , k , and k , correspond to a half-life of 31 and 62 hours respectively and are chosen toachieve a realistic model outcome. Consequently generation rate constants k , and k , arechosen such that a surface traction of 0 . .
25 kN over asurface of 300mm – results in an equilibrium state for a = a = 1 when c = 1, that is whenthe concentration of osteoblast equals that of healthy bone.(e) The diffusivity D ( ρ ) = k (1 − ρ ) is controlled by the porosity 1 − ρ and the constant k . With k = 260mm / month we set it to a standard value for the diffusion of bioactive molecules thatis measured for soluble proteins, see [5] and [51].(f) We use Voigt’s bound as an approximation of the material properties of the bone-scaffoldcomposite. More precisely, we model bone and PCL as linear isotropic materials with material N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 17
Figure 1.
Experiment including fixateur. Shown is a vertical section through thecylindrical defect site. Fixateur domain is colored in gold. From left to right: regen-erated bone at 3 months, 12 months and a view on top of the defect site. The greycolored areas illustrate the top and bottom cylinder/fixateur caps.constants chosen as collected in Table 1. The effective properties of the compositum are thenobtained by adding the weighted tensors.(g) The constant k drives the rate of bone regeneration, k is related to the overall osteoblastproduction and k influences the effect of osteoblast proliferation. These values are fitted toachieve realistic outcome in the simulations. Table 1.
Parameters for the bone regeneration model
Param. Value Description T
12 months Period of bone regenerationΩ L = 30mm, r = 10mm Cylinder with length L , radius rρ ρ ≡ .
13 Scaffold volume fraction C ( ρ, σ, b ) b C b + ρσ C ρ Voigt bound for composites D ( ρ ) k (1 − ρ ) Diffusivity of bioactive molecules( λ b , µ b ) (2.88GPa, 1.92GPa) Derived from ( E b , ν b ) = (5 GPa , . λ ρ , µ ρ ) (1.97GPa, 0.17GPa) Derived from ( E ρ , ν ρ ) = (0 . , . C b C b A = λ b tr( A ) Id +2 µ b A Material tensor of healthy bone C ρ C ρ A = λ ρ tr( A ) Id +2 µ ρ A Material tensor of PCL k . k , k , k ,
16 Decay rate first molecule k , k . k
260 mm / month Diffusivity of the a i w/o scaffold k . k .
07 Proliferation constant for osteoblasts4.2.
Numerical Implementations.
We use a simple first-order implicit in time Euler scheme to solvethe equations displayed in the order displayed above. The fact that an implicit approach is feasibleis due to the simple structure of the ODEs and the linearity of the diffusion equation. It is worth
Figure 2.
Experiment excludig fixateur. Shown is a vertical section through thecylindrical defect site. From left to right: regenerated bone at 3 months, 12 monthsand a view on top of the defect site. The grey colored areas illustrate the top andbottom cylinder caps.mentioning that this reduces the computational cost of solving the system drastically as only very fewtime steps are needed to achieve acceptable accuracy in the simulations. The elastic and the diffusionequation are discretized using P1 elements and the meshes were generated using the ComputationalGeometry Algorithms Library CGAL [6].4.3.
Discussion of Numerical Simulations.
In Figure 1 the domain of computation with an addedfixateur is shown. Here we assume the material of the fixateur to be titanium with Young’s moduluschosen to 100GPa and a Poisson’s ratio of 0.31. Bone growth and osteoblast production is disabled x (mm) b b after 12 months b after 3 months Figure 3.
Relative bone density averagedover horizontal slices after 3 and 12 monthsin the experiment including the fixateur. in the space occupied by the fixateur. In Figure3 we present the relative bone density averagedover horizontal slices in the fixateur experimentat 3 and 12 months. We observe that both the re-generated bone after 3 and after 12 months agreewell with the experimental results shown in [15,Figure 2, ‘Scaffold only’]. There, the same shapeof regenerated bone, with a flat area in the mid-dle of the defect site and a significant gradienttowards the proximal and distal interface, is ob-served.In Figure 1, the result of the stress shielding ef-fect of the fixateur is clearly visible, with littleregenerated bone in the central part of the defectsite close to the fixateur. Comparing to [49, Figs4C, 5C] or [41, Figures 3a, 3b] we see that thisis also observed in experiments. Bone mass lossdue to stress shielding is indeed a long recognized,major issue in orthopaedic surgery [45, 48].The computation excluding the fixateur is per-formed using a reduced surface traction that isset to 70% of the surface traction in the fixateurmodel to account for the stress shielding of the fixateur. This experiment is the direct analogon ofthe 1D model in [39]. Naturally, we see that bone regenerates symmetrically and that the result is
N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 19 essentially a one dimensional distribution of bone comparable to the results in [39]. Note that theasymmetries encountered in the more realistic model including the fixateur can not be resolved bya one-dimensional simplification. This has important implications for the porosity optimization ofscaffolds where a three dimensional simulation can thus help to achieve a more appropriate optimaldesign.
Acknowledgements
The authors gratefully acknowledge support from BMBF within the e:Med program in the SyMBoDconsortium (grant number 01ZX1910C). Furthermore the authors thank Luca Courte (Freiburg) andDorothee Knees (Kassel) for helpful suggestions and discussions.
Appendix A. Properties of Diffusion Equation
This section provides the regularity results needed in the existence proof of Theorem 3.2. We begin bystating the regularity results in section A.1, then discuss the notion of Gr¨oger regular sets in sectionA.2 and conclude with the proofs in section A.3.A.1.
Regularity Results.
Let Ω be a Lipschitz domain with a Dirichlet-Neumann partition ofthe boundary ∂ Ω = Γ N ∪ Γ D . Both Γ D and Γ N are allowed to have vanishing measure. Let D ∈ L ∞ (Ω , R n × n ) be uniformly elliptic, k > f , f ∈ L ( I, L (Ω)) and a ∈ L ∞ (Ω), then we are in-terested in L ( I, C (Ω)), L ∞ ( I, L ∞ (Ω)) and L ( I, C α (Ω)) regularity of a i ∈ H ( I, H D (Ω) , H D (Ω) (cid:48) ), i = 1 , (cid:90) I (cid:104) d t a i , · (cid:105) H D (Ω) dt + (cid:90) I (cid:90) Ω D ∇ a ∇ · + ka · dxdt = (cid:90) I (cid:90) Ω f i · dxdt (A.1) a i (0) = a . (A.2)We are also interested in the regularity of the difference a − a . Note that a − a has better regularityproperties as the initial value a (0) − a (0) = 0 is smooth. We will need varying assumptions in additionto the ones stated above, depending on the regularity we are after. Note that the main difficulty stemsfrom the mixed boundary conditions as in this case the usual elliptic H (Ω) regularity fails, see forexample [44, 26, 21]. Let us now state our main theorems. Theorem A.1 ( L ( I, C α (Ω)) Regularity) . Assume that Ω ⊂ R n with n = 1 , , is a Lipschitz domain, Γ N ∪ Γ D = ∂ Ω and assume that Ω ∪ Γ N is Gr¨oger regular, see A.4. Furthermore, let f i ∈ L ( I, L (Ω)) , D ∈ L ∞ (Ω , M s ) , k > and a ∈ L ∞ (Ω) . Then there is α > such that it holds a − a ∈ L ( I, C α (Ω)) and (A.3) (cid:107) a − a (cid:107) L ( I,C α (Ω)) ≤ C (cid:107) f − f (cid:107) L ( I,L (Ω)) , In addition, for every ε > it holds a i ∈ L ([ ε, T ] , C α (Ω)) and a i ∈ L ( I, C (Ω)) . Theorem A.2.
Assume that f ∈ L ∞ ( I × Ω) , a ∈ L ∞ (Ω) and the assumptions of the beginning ofthe section. Then it holds that a i ∈ L ∞ ( I × Ω) . The above regularity theorems apply to ˜ a i of the main body of the article, where ˜ a i is the part of thesolution corresponding to homogeneous boundary conditions. Clearly all the results still hold true for˜ a i + 1. To conclude we need positivity of ˜ a i + 1, therefore we consider a slightly different equation than(A.1). Theorem A.3 (Positivity) . Assume a function a ∈ H ( I, H D (Ω) , H D (Ω) (cid:48) ) satisfies the followingequality in the space L ( I, H D (Ω)) (cid:48) (cid:90) I (cid:104) d t a, · (cid:105) H D (Ω) dt + (cid:90) I (cid:90) Ω D ∇ a ∇ · + k ( a + 1) · dxdt = (cid:90) I (cid:90) Ω f · dxdt, and a (0) + 1 = 0 . Furthermore suppose that f ∈ L ( I, L (Ω)) is non-negative and the remainingassumptions stated at the beginning of the appendix hold true. Then a + 1 ≥ . A.2.
Boundary Regularity.
We say a bounded, open set Ω ⊂ R n is a Lipschitz domain if Ω is aLipschitz manifold with boundary, see [21, Definition 1.2.1.2]. In the following we will denote the cube[ − , n ⊂ R n by Q . The following definition is due to Gr¨oger, see [22]. Definition A.4 (Gr¨oger Regular Sets) . Let Ω ⊂ R n be bounded and open and Γ ⊂ ∂ Ω a relativelyopen set. We call Ω ∪ Γ Gr¨oger regular, if for every x ∈ ∂ Ω there are open sets
U, V ⊂ R n with x ∈ U , and a bijective, bi-Lipschitz map φ : U → V , such that φ ( x ) = 0 and φ ( U ∩ (Ω ∪ Γ)) is either { x ∈ Q | x n < } , { x ∈ Q | x n ≤ } or { x ∈ Q | x n < } ∪ { x ∈ Q | x n = 0 , x n − < } .It can easily be seen that a Gr¨oger regular set Ω (no matter the choice Γ ⊂ ∂ Ω) is a Lipschitzdomain, see [24, Theorem 5.1]. The next two theorems characterize Gr¨oger regular sets in two andthree dimension. We cite the results from [24].
Theorem A.5 (Gr¨oger Regular Sets in 2D) . Let Ω ⊂ R be a Lipschitz domain and Γ ⊂ ∂ Ω berelatively open. Then Ω ∪ Γ is Gr¨oger regular if and only if Γ ∩ ( ∂ Ω \ Γ) is finite and no connectedcomponent of ∂ Ω \ Γ consists of a single point. Theorem A.6 (Gr¨oger Regular Sets in 3D) . Let Ω ⊂ R be a Lipschitz domain and Γ ⊂ ∂ Ω berelatively open. Then Ω ∪ Γ is Gr¨oger regular if and only if the following two conditions hold(i) ∂ Ω \ Γ is the closure of its interior.(ii) For any x ∈ Γ ∩ ( ∂ Ω \ Γ) there is an open neighborhood U x of x and a bi-Lipschitz map φ : U x ∩ Γ ∩ ( ∂ Ω \ Γ) → ( − , . A.3.
Proofs of the Regularity Results.
We will begin with the L ( I, C α (Ω)) results, i.e., the proofof theorem A.1. To this end we need two results from the literature, one, [24, Theorem 3.3] on mixedelliptic boundary value problems that will yield the C α (Ω) information as well as a maximal L (Ω)regularity result from [4] to transfer this to the solution of the associated parabolic equation. Theorem A.7 (H¨older Regularity) . Let Ω ⊂ R n for n = 2 , , be open, bounded and connected andassume that Ω ∪ Γ N is Gr¨oger regular. Assume further, that D ∈ L ∞ (Ω , R n × n ) is uniformly elliptic, k > and let q > n and denote by q (cid:48) its adjoint exponent, i.e., /q + 1 /q (cid:48) = 1 . Define the operator M : H D (Ω) → H D (Ω) (cid:48) with M v = (cid:90) Ω D ∇ v ∇ · + kv · dx. Then there is α > such that M − : W ,q (cid:48) D (Ω) (cid:48) → C α (Ω) is well defined and continuous.Proof. See [24, Theorem 3.3]. (cid:3)
In particular, in dimensions n = 2 , q = 4 and by the Sobolev embedding theorems,see [21], it holds that L (Ω) (cid:44) → W ,q (cid:48) D (Ω) (cid:48) and hence for f ∈ L (Ω) a solution v ∈ H D (Ω) to (cid:90) Ω D ∇ v ∇ · + kv · dx = (cid:90) Ω f · dx in H D (Ω) (cid:48) lies in C α (Ω) and satisfies the estimate (cid:107) v (cid:107) C α (Ω) ≤ C (cid:107) f (cid:107) L (Ω) . In order to amplify this elliptic regularity result to the parabolic setting, we will use an L (Ω) maximalregularity result. In our particular case this is a statement of the form: f ∈ L ( I, L (Ω)) implies thatthe solution a ∈ H ( I, H D (Ω) , H D (Ω) (cid:48) ) of d t a + M a = f and a (0) = a is of regularity H ( I, L (Ω)). Consequently f , d t a and M a are members of the space L ( I, L (Ω)),hence the term maximal regularity. Depending on the operator M , problems of this form can be N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 21 delicate, we refer to [31] for certain results and to [4] for a recent discussion. In our case M is self-adjoint and does not depend on time and therefore L (Ω) maximal regularity holds Theorem A.8 (Maximal L (Ω) Regularity) . Let Ω ⊂ R n be a Lipschitz domain, ∂ Ω = Γ N ∪ Γ D , D ∈ L ∞ (Ω , M s ) uniformly elliptic and symmetric and k > . Define M : H D (Ω) → H D (Ω) (cid:48) with M v = (cid:90) Ω D ∇ v ∇ · + kv · dx. and also its part in L (Ω) which we will denote by M , i.e., dom( M ) := { v ∈ H D (Ω) | M v ∈ L (Ω) } with M v = M v Then set X := H ( I, L (Ω)) ∩ L ( I, dom( M )) endowed with the norm (cid:107) a (cid:107) X := (cid:107) d t a (cid:107) L ( I,L (Ω)) + (cid:107) M a ( · ) (cid:107) L ( I,L (Ω)) . Here H ( I, L (Ω)) shall denote the functions in H ( I, L (Ω)) with vanishing initial value. Then thespace X is a Hilbert space and the map d t + M : X → L ( I, L (Ω)) with a (cid:55)→ d t a + M a is a linear homeomorphism.
Addendum.
Furthermore, if we consider the problem as above but arbitrary initial value in H D (Ω) ,i.e., d t a + M a = f and a (0) = a ∈ H D (Ω) , the solution a lies in H ( I, L (Ω)) .Proof. This is discussed in [4, Section 4]. The Addendum requires a to lie in trace space (the spaceof initial values) of H ( I, L (Ω)) ∩ L ( I, dom( M )). As we assume D to be symmetric we have that M is self-adjoint and hence H D (Ω) = dom( M / ) which coincides with the trace space, see [4], sections 2and 4. (cid:3) Proof of Theorem A.1.
We will first show that a − a ∈ L ( I, C α (Ω)) and that the estimate (A.3) issatisfied. Note that a := a − a satisfies an equation of type (A.1) with right hand side f := f − f ∈ L ( I, L (Ω)) and zero initial condition, namely d t a + M a = f in L ( I, H D (Ω)) (cid:48) and a (0) = 0 . by Theorem A.8 a lies in the space X and satisfies almost everywhere in I the equation d t a ( t )+ M a ( t ) = f ( t ). Using Theorem A.7 we can estimate(A.4) (cid:107) a ( t ) (cid:107) C α (Ω) ≤ C (cid:107) f ( t ) − d t a ( t ) (cid:107) L (Ω) . It follows (cid:107) a (cid:107) L ( I,C α (Ω)) ≤ C ( (cid:107) f (cid:107) L ( I,L (Ω)) + (cid:107) d t (cid:107) L ( X ,L ( I,L (Ω))) (cid:107) a (cid:107) X )and using the continuity of ( d t + M ) − established in Theorem A.8 we get the estimate (cid:107) a (cid:107) X ≤ C (cid:107) f (cid:107) L ( I,L (Ω)) which yields combined with the previous estimates (cid:107) a − a (cid:107) L ( I,C α (Ω)) ≤ C (cid:107) f (cid:107) L ( I,L (Ω)) . Now we show that a i is a member of L ( I, C (Ω)). We decompose a i into a i = a i + a i where a i solves d t a i + M a i = f i and a i (0) = 0and a i solves d t a i + M a i = 0 and a i (0) = a . By repeating our previous computations it is clear that a i ∈ L ( I, C α (Ω)). To conclude we will provethat for every ε > a i ∈ L ([ ε, T ] , C α (Ω)) ∩ L ∞ ( I × Ω) . Let us begin with a i ∈ L ([ ε, t ] , C α (Ω)). As a i ∈ C ( I, L (Ω)) ∩ L ( I, H D (Ω)) point evaluations arewell defined and there is a sequence ( ε n ) of real numbers with ε n → a i ( ε n ) ∈ H D (Ω) forall n ∈ N . For given ε > ε n such that ε n ≤ ε . Then apply the Addendum of Theorem A.8 toobtain a i ∈ H ([ ε n , T ] , L (Ω)) ∩ L ([ ε n , T ] , dom( M )) . Repeating the computations in the beginning of this proof we arrive at (A.4) with a i instead of a . Itis then clear that a i ∈ L ([ ε, T ] , C α (Ω)). The L ∞ ( I × Ω) part of the theorem holds more generally forright hand sides in L ∞ ( I × Ω) and is addressed in the following proof of Theorem A.2. (cid:3)
Proof of Theorem A.2.
We will use Stampacchias truncation method [46], that is for a real number ¯ a we will test the PDE with( a i − ¯ a ) + := max(0 , a i − ¯ a ) and ( a i + ¯ a ) − := − min(0 , a i + ¯ a ) . One can show that ( a i − ¯ a ) + and ( a i + ¯ a ) − are members of H ( I, H D (Ω) , H D (Ω) (cid:48) ) if a i is itself in thatspace and that it holds (cid:90) t (cid:104) d t a i , ( a i − ¯ a ) + (cid:105) H D (Ω) dt = 12 (cid:13)(cid:13) ( a i − ¯ a ) + ( t ) (cid:13)(cid:13) L (Ω) − (cid:13)(cid:13) ( a − ¯ a ) + (cid:13)(cid:13) L (Ω) and (cid:90) Ω D ∇ a i ( t ) ∇ ( a i − ¯ a ) + ( t ) dx = (cid:90) Ω D ∇ ( a i − ¯ a ) + ( t ) ∇ ( a i − ¯ a ) + ( t ) dx such as (cid:90) Ω a i ( t )( a i − ¯ a ) + ( t ) dx = (cid:13)(cid:13) ( a i − ¯ a ) + ( t ) (cid:13)(cid:13) L (Ω) + (cid:90) Ω ¯ a ( a i − ¯ a ) + ( t ) dx Hence it follows for every t ∈ I and ¯ a ≥ max( (cid:107) f (cid:107) L ∞ ( I × Ω) , (cid:107) a (cid:107) L ∞ (Ω) )12 (cid:13)(cid:13) ( a i − ¯ a ) + ( t ) (cid:13)(cid:13) L (Ω) ≤ (cid:90) t (cid:90) Ω ( f − ¯ a )( a i − ¯ a ) + dxds + 12 (cid:13)(cid:13) ( a − ¯ a ) + (cid:13)(cid:13) L (Ω) ≤ . This implies that a i ≤ ¯ a almost everywhere in I × Ω. Similarly one establishes ¯ a ≤ a i using the testfunction ( a i + ¯ a ) − . (cid:3) Proof of Theorem A.3.
One can check that ( a +1) − = − min(0 , a +1) is a member of H ( I, H D (Ω) , H D (Ω) (cid:48) )and that it holds for all t ∈ I (cid:90) t (cid:104) d t a, ( a + 1) − (cid:105) H D ds = 12 (cid:16) (cid:13)(cid:13) ( a + 1) − (0) (cid:13)(cid:13) L (Ω) − (cid:13)(cid:13) ( a + 1) − ( t ) (cid:13)(cid:13) L (Ω) (cid:17) and (cid:90) t (cid:90) Ω (cid:104) D ∇ a, ∇ ( a + 1) − (cid:105) + k ( a + 1)( a + 1) − dxds = − (cid:90) t (cid:90) Ω (cid:104) D ∇ (cid:2) ( a + 1) − (cid:3) , ∇ (cid:2) ( a + 1) − (cid:3) (cid:105) + k (cid:2) ( a + 1) − (cid:3) dxds ≤ . Testing the full equation with ( a + 1) − and using these computations one finds12 (cid:13)(cid:13) ( a + 1) − ( t ) (cid:13)(cid:13) L (Ω) ≤ − (cid:90) t (cid:90) Ω f · ( a + 1) − dxds ≤ . (cid:3) N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 23
Appendix B. Ordinary Differential Equations
The plan for this section is as follows. We recall Gr¨onwall’s inequality in Theorem B.1 and will thereafterprovide a Banach space valued ODE existence theorem in Theorem B.2. The version of our ODEtheorem guarantees existence and uniqueness of short-time solutions. We also show that the solutionsdepend continuously on the initial value with a common short-time existence interval, at least locallyaround a fixed initial value. This will come in handy for analyzing pointwise properties of certain ODEsin Lemma B.5. Finally we prove the existence and uniqueness of the cell ODE in Lemma B.6 and thebone ODE in (B.7)
Lemma B.1 (Gr¨onwall Variant) . Let X be a Banach space and I = [0 , T ] an interval. Let x , y ∈ X and assume that γ , γ are members of L ( I, X ) . Define x and y to be the integral curves x ( t ) = x + (cid:90) t γ ( s ) ds and y ( t ) = y + (cid:90) t γ ( s ) ds. Now assume that we can estimate the integrants γ , γ in the following form (cid:107) γ ( t ) − γ ( t ) (cid:107) ≤ α ( t ) + β ( t ) (cid:107) x ( t ) − y ( t ) (cid:107) for all t ∈ I, where α, β ∈ L ( I ) are non-negative functions. Then it holds that (cid:107) x ( t ) − y ( t ) (cid:107) ≤ C (cid:0) (cid:107) x − y (cid:107) + (cid:107) α (cid:107) L ( I ) (cid:1) for all t ∈ I and the constant C can be chosen to be C = 1 + (cid:107) β (cid:107) L ( I ) exp( (cid:107) β (cid:107) L ( I ) ) .Proof. Just write out the estimate that the difference (cid:107) x ( t ) − y ( t ) (cid:107) satisfies due to the assumptionsand then use the usual integral formulation of Gr¨onwall’s inequality, see for example [40, Theorem1.2.8]. (cid:3) Theorem B.2 (Local Existence) . Let X be a Banach space, I = [ a, b ] a bounded interval and F : I × X → X a Carath´eodory function, i.e., F ( · , x ) is Bochner measurable for all x ∈ X and F ( t, · ) is continuous almost everywhere in I . Assume that for every bounded set B ⊂ X there are functions m B ∈ L p ( I ) , p ∈ [1 , ∞ ] and L B ∈ L ( I ) , possibly depending on B , such that (cid:107) F ( t, x ) (cid:107) X ≤ m B ( t ) a.e. in I, ∀ x ∈ B, (B.1) (cid:107) F ( t, x ) − F ( t, y ) (cid:107) X ≤ L B ( t ) (cid:107) x − y (cid:107) a.e. in I, ∀ x, y ∈ B. (B.2) Let furthermore t ∈ I , x ∈ X and R > be arbitrary, then there exists a time interval I δ :=[ t − δ, t + δ ] ∩ I such that for any initial value y ∈ B R ( x ) ⊂ X there is a unique short time solution x ∈ W ,p ( I δ , X ) of the ODE d t x ( t ) = F ( t, x ( t )) and x ( t ) = y . Moreover the map y (cid:55)→ x ( y ) taking the initial value to its solution seen as a map B R ( x ) ⊂ X → C ( I δ , X ) is continuous.Remark B.3 . Note that the Lipschitz assumption (B.2) implies that x (cid:55)→ F ( t, x ) is continuous. There-fore, to establish the Carath´eodory regularity of F it is enough to provide the Bochner measurabilityof the maps F ( · , x ) : I → X for all x ∈ X. Proof.
First we clarify the dependence of δ ; we choose it to satisfy (cid:90) I δ m B R ( x ) ( s ) ds ≤ R and (cid:90) I δ L B R ( x ) ( s ) ds < . (B.3)Then we consider the complete metric space E ⊂ C ( I δ , X ) given by E := { x ∈ C ( I δ , X ) | sup t ∈ I δ (cid:107) x ( t ) − x (cid:107) ≤ R } and define the map Φ : E → E with Φ( x )( t ) = y + (cid:90) tt F ( s, x ( s )) ds. We now proceed by showing the following facts(i) For all x ∈ E the map t (cid:55)→ F ( t, x ( t )) is Bochner integrable as a map I δ → X . In fact it is amember of L p ( I δ , X ).(ii) The function Φ is a self-mapping and a contraction.(iii) The fix-point of Φ is a member of W ,p ( I δ , X ) and corresponds to the solution of the ODE.(iv) The solution depends continuously on the initial data.To establish (i) note that the assumption of Carath´eodory regularity of F implies that t (cid:55)→ F ( t, x ( t ))is Bochner measurable as a map I δ → X for all maps x : I δ → X that are itself Bochner measurable,which clearly holds for members of E . We are left to show the integrability, so we estimate for x ∈ E (cid:90) I δ (cid:107) F ( t, x ( t )) (cid:107) pX dt ≤ (cid:90) I δ | m B R ( x ) ( t ) | p dt < ∞ , which shows the assertion. Now to (ii). Let again x ∈ E and estimatesup t ∈ I δ (cid:107) Φ( x )( t ) − x (cid:107) X ≤ (cid:107) x − y (cid:107) + (cid:90) I δ (cid:107) F ( t, x ( t )) (cid:107) X dt ≤ R + (cid:90) I δ m B R ( x ) ( t ) dt ≤ R. To see that Φ is a contraction compute for x, y ∈ E sup t ∈ I δ (cid:107) Φ( x )( t ) − Φ( y )( t ) (cid:107) X ≤ sup t ∈ I δ (cid:90) I δ (cid:107) F ( t, x ( t )) − F ( t, y ( t )) (cid:107) X dt ≤ (cid:13)(cid:13) L B R ( x ) (cid:13)(cid:13) L ( I δ ) (cid:124) (cid:123)(cid:122) (cid:125) < (cid:107) x − y (cid:107) E The claim (iii) follows as the unique fix-point x of Φ is, by the fundamental theorem, a solution to theODE. As d t x ( t ) = F ( t, x ( t )) the L p integrability of the derivative of this fix-point follows from the oneof F ( · , x ( · )) which was established in (i).Finally to (iv), where we will employ Gr¨onwall’s lemma. Let y and y be in B R ( x ), then the accordingsolutions are given by y ( t ) = y + (cid:90) tt F ( s, y ( s )) ds and y ( t ) = y + (cid:90) tt F ( s, y ( s )) ds. The difference in the integrands can be estimated by (cid:107) F ( t, y ( t )) − F ( t, y ( t )) (cid:107) X ≤ CL B R ( x ) ( t ) (cid:107) y ( t ) − y ( t ) (cid:107) X . So applying Lemma B.1 with α = 0 and β = L B R ( x ) yields (cid:107) y ( t ) − y ( t ) (cid:107) X ≤ C (cid:107) y − y (cid:107) X . (cid:3) Remark
B.4 . It is often of interest to show the existence of long-time solutions. A particularly simplecase in the setting of the above theorem is encountered if it holds that for any initial value y ∈ B R ( x ) ⊂ X the solution takes values only in B R ( x ). Then one glues together multiple short-timesolutions with the guarantee of δ not deteriorating. N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 25
Lemma B.5 (Pointwise Properties of Solutions) . Let I = [ a, b ] be an interval and K : I × R → R aCarath´eodory function such that x ≥ implies K ( t, x ) ≥ for all t ∈ I . For fixed t ∈ I consider theODE x (cid:48) ( t ) = K ( t, x ( t )) (cid:18) − x ( t ) θ (cid:19) and x ( t ) = x , where θ > and λ ≥ are fixed numbers and x ∈ [0 , θ ] . Assume that there is an interval I δ =[ t − δ, t + δ ] ∩ I such that for any choice of x ∈ [0 , θ ] we have a solution x ∈ W , ( I δ ) of the ODEwhich we assume to continuously depend on the initial data x , i.e., we assume that for every x ∈ [0 , θ ] there is a neighborhood N x around x such that x (cid:55)→ x is continuous as a map N x → C ( I δ ) , where x is the solution to the ODE with initial value x . Then it holds ≤ x ( t ) ≤ θ for all t ∈ I δ . Proof.
We know that x satisfies the identity x ( t ) = x + (cid:90) tt K ( s, x ( s )) (cid:18) − x ( s ) θ (cid:19) ds for all t ∈ I δ Upper Barrier.
We prove this by contradiction. Suppose there was s ∈ I δ with x ( s ) > θ , then on aneighborhood of s solution must be non-increasing which can be seen as follows: Due to the continuityof x there is ε > x ( t ) ≥ θ for all t ∈ ( s − ε, s + ε ) . If x was not non-increasing on ( s − ε, s + ε ) then there exist t < t in that interval such that x ( t ) > x ( t ) and therefore0 < x ( t ) − x ( t ) = (cid:90) t t K ( s, x ( s )) (cid:124) (cid:123)(cid:122) (cid:125) ≥ (cid:18) − x ( s ) θ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) ≤ ds ≤ , which settles the claim. Now, by judicious Zornification we produce a maximal interval Z around s on which x is non-increasing. Then t ∗ := inf Z = t (if it was not t , repeat the above reasoning andfind that Z was not maximal) and hence θ < x ( t ∗ ) = x ( t ) ≤ θ clearly is a contradiction. Lower Barrier.
With an analogue reasoning as in the proof for the upper barrier we can establishthe following: If x ( s ) ∈ (0 , θ ] for some s ∈ I δ then x ( t ) ∈ [ x ( s ) , θ ] for all t ≥ s . This yields the claimfor all initial values strictly larger than zero. We need x ( s ) to exceed zero to guarantee the existenceof a small interval ( s − δ, s + δ ) where x ≥ K ( s, x ( s )) ≥ x ( t ) = 0 we approximate the solution by considering initial values x n ( t ) = 1 /n , i.e., wefind solutions x n to x (cid:48) n ( t ) = K ( s, x n ( s )) (cid:18) − x n ( s ) θ (cid:19) with x n ( t ) = n − . As shown before we then know that x n ( t ) ∈ [1 /n, θ ] for all t ∈ I δ . By the continuity we assumed wecan pass to the limit in n and obtain 0 ≤ x ( t ) ≤ θ for all t ∈ I δ . (cid:3) We continue by explaining the connection between Banach space valued ODEs and the formulationas a family of real valued ODEs in our examples (2.13), (2.15). A moments reflection reveals that weonly need to guarantee that for every fixed x ∈ Ω it holds t (cid:55)→ d t b ( t, x ) = ddt ( t (cid:55)→ b ( t, x )) , where on the right hand side we talk about the usual, real valued, weak derivative. Therefore let b ∈ W , ( I, C (Ω)), then for d t b it holds by definition (cid:90) I d t b ( t ) ϕ ( t ) dt = − (cid:90) I b ( t ) ∂ t ϕ ( t ) dt ∀ ϕ ∈ D ( I ) . The integral used above is the C (Ω) valued Bochner integral, thus using that point evaluation is alinear and continuous map on the space of continuous functions we find that for every x ∈ Ω it holds (cid:90) I d t b ( t )( x ) ϕ ( t ) dt = − (cid:90) I b ( t )( x ) ∂ t ϕ ( t ) dt ∀ ϕ ∈ D ( I ) , meaning that for every fixed x ∈ Ω the function t (cid:55)→ b ( t )( x ) satisfies the real valued ODE as desired.If we additionally assume that both, the Banach space valued ODE and the scalar one have uniquesolutions, we obtain that both settings are equivalent. This means that the above Lemma is applicableto deduce pointwise properties of the solutions to the Banach space valued ODEs. Lemma B.6 (Solveability of the Cell ODE) . Assume H : R N +2 → R is locally Lipschitz continuousand that H is non-negative if all its arguments are non-negative. Further, let the assumptions 3.1 besatisfied. Then there is a unique solution c ∈ W ,p ( I, C (Ω)) of the equation d t c = H ( a , . . . , a N , b, c ) (cid:18) − c − ρ (cid:19) with c (0) = 0 . Furthermore the solution satisfies ≤ c ( t, x ) ≤ − ρ ( x ) for all t ∈ I and x ∈ Ω .Proof. To begin with, define the auxiliary function˜ H : R N +2 × [ c P , C P ] → R with ˜ H ( a, b, c, ρ ) = H ( a, b, c ) (cid:18) − c − ρ (cid:19) . Then ˜ H is locally Lipschitz continuous. Now note that the ODE is induced by F : I × C (Ω) → C (Ω) with F ( t, c ) = x (cid:55)→ ˜ H ( a ( t, x ) , b ( x ) , c ( x ) , ρ ( x )) . We aim to apply Theorem B.2 to produce a short-time solution, hence we need to guarantee(i) F ( t, c ) ∈ C (Ω) for all t ∈ I and c ∈ C (Ω),(ii) F ( · , c ) : I → C (Ω) is Bochner measurable for all c ∈ C (Ω),(iii) F satisfies (B.1) and (B.2).The statement (i) is clear as F ( t, c ) is a composition of continuous functions. To prove (ii) we write a asa pointwise almost everywhere limit of finitely valued, measurable functions ( s k ) ⊂ S ( I, C (Ω) N ). Then t (cid:55)→ ˜ H ( s k ( t ) , b, c, ρ ) is still a member of S ( I, C (Ω) N ). As s k ( t ) → a ( t ) in C (Ω) almost everywhere in I , for fixed t ∈ I the set (cid:91) k ∈ N { ( s k ( t, x ) , b ( x ) , c ( x ) , ρ ( x )) , ( a ( t, x ) , b ( x ) , c ( x ) , ρ ( x )) | x ∈ Ω } ⊂ R N +3 is relatively compact in R N +3 . Hence it holds (cid:13)(cid:13)(cid:13) ˜ H ( s k ( t ) , b, c, ρ ) − ˜ H ( a ( t ) , b, c, ρ ) (cid:13)(cid:13)(cid:13) C (Ω) ≤ C (cid:107) s k ( t ) − a ( t ) (cid:107) C (Ω) . This establishes the Bochner measurability in ( ii ). To show (iii) let B ⊂ C (Ω) be bounded andcompute (cid:107) F ( t, ˜ c ) (cid:107) C (Ω) ≤ (cid:107) H ( a ( t ) , b ( t ) , ˜ c ) (cid:107) C (Ω) (cid:13)(cid:13)(cid:13)(cid:13) − ˜ c − ρ (cid:13)(cid:13)(cid:13)(cid:13) C (Ω) ≤ m HB ( t ) (cid:0) x ∈ Ω [1 / (1 − ρ ( x ))] (cid:107) ˜ c (cid:107) C (Ω) (cid:1) ≤ Cm HB ( t ) . N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 27
In a similar way we estimate (cid:13)(cid:13) F ( t, ˜ c ) − F ( t, ˜˜ c ) (cid:13)(cid:13) C (Ω) ≤ (cid:13)(cid:13) H ( a ( t ) , b ( t ) , ˜ c ) − H ( a ( t ) , b ( t ) , ˜˜ c ) (cid:13)(cid:13) C (Ω) (cid:13)(cid:13)(cid:13)(cid:13) − ˜ c − ρ (cid:13)(cid:13)(cid:13)(cid:13) + (cid:107) H ( a ( t ) , b ( t ) , ˜ c ) (cid:107) C (Ω) (cid:13)(cid:13)(cid:13)(cid:13) ˜ c − ρ − ˜˜ c − ρ (cid:13)(cid:13)(cid:13)(cid:13) C (Ω) ≤ CL HB ( t ) (cid:13)(cid:13) ˜ c − ˜˜ c (cid:13)(cid:13) C (Ω) + Cm HB ( t ) (cid:13)(cid:13) ˜ c − ˜˜ c (cid:13)(cid:13) C (Ω) ≤ C max( L HB ( t ) , m HB ( t )) (cid:124) (cid:123)(cid:122) (cid:125) = ζ (cid:13)(cid:13) ˜ c − ˜˜ c (cid:13)(cid:13) C (Ω) and ζ is a member of L ( I ). To establish a long-time solution note that by our pointwise lemma B.5we have 0 ≤ c ( t, x ) ≤ − ρ ( x ) ≤ . Remember that we discussed the connection between Banach space valued ODEs and R valued ODEsin the section 2. Finally using the remark following Theorem B.2 we conclude. (cid:3) Clearly the bone ODE can now be treated identically, provided one assumes the same for K as onedid for H . Lemma B.7 (Solveability of the Bone ODE) . Assume K : R N +2 → R is locally Lipschitz continuousand that K is non-negative if all its arguments are non-negative. Further, let the assumptions 3.1 besatisfied. Then there is a unique solution b ∈ W ,q ( I, C (Ω)) of the equation d t b = H ( a , . . . , a N , b, c ) (cid:18) − b − ρ (cid:19) with b (0) = 0 . Furthermore the solution satisfies ≤ b ( t, x ) ≤ − ρ ( x ) for all t ∈ I and x ∈ Ω . References [1] R. A. Adams and J. J. Fournier.
Sobolev spaces . Elsevier, 2003.[2] J. Alierta, M. P´erez, and J. Garc´ıa-Aznar. An interface finite element model can be used to predict healing outcomeof bone fractures.
Journal of the Mechanical Behavior of Biomedical Materials , 29:328 – 338, 2014.[3] G. Allaire.
Shape optimization by the homogenization method , volume 146. Springer Science & Business Media,2012.[4] W. Arendt, D. Dier, and S. Fackler. JL Lions’ problem on maximal regularity.
Archiv der Mathematik , 109(1):59–72,2017.[5] A. Badugu, C. Kraemer, P. Germann, D. Menshykau, and D. Iber. Digit patterning during limb development as aresult of the BMP-receptor interaction.
Scientific reports , 2:991, 2012.[6] J.-D. Boissonnat, O. Devillers, M. Teillaud, and M. Yvinec. Triangulations in CGAL. In
Proceedings of the sixteenthannual symposium on Computational geometry , pages 11–18, 2000.[7] F. Boyer and P. Fabrie.
Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and RelatedModels , volume 183. Springer Science & Business Media, 2012.[8] H. Brezis.
Functional analysis, Sobolev spaces and partial differential equations . Springer Science & Business Media,2010.[9] G. M. Calori, E. L. Mazza, S. Mazzola, A. Colombo, F. Giardina, F. Roman`o, and M. Colombo. Non-unions.
ClinicalCases in Mineral and Bone Metabolism , 14(2):186, 2017.[10] V. J. Challis, J. K. Guest, J. F. Grotowski, and A. P. Roberts. Computationally generated cross-property boundsfor stiffness and fluid permeability using topology optimization.
International Journal of Solids and Structures ,49(23-24):3397–3408, 2012.[11] S. Checa and P. J. Prendergast. Effect of cell seeding and mechanical loading on vascularization and tissue for-mation inside a scaffold: A mechano-biological model using a lattice approach to simulate cell activity.
Journal ofBiomechanics , 43(5):961 – 968, 2010.[12] P. G. Ciarlet.
Mathematical Elasticity: Volume I: three-dimensional elasticity . North-Holland, 1988.[13] P. G. Ciarlet. On Korn’s inequality.
Chinese Annals of Mathematics, Series B , 31(5):607–618, 2010.[14] A. Cipitria, C. Lange, H. Schell, W. Wagermaier, J. C. Reichert, D. W. Hutmacher, P. Fratzl, and G. N. Duda.Porous scaffold architecture guides tissue formation.
Journal of Bone and Mineral Research , 27(6):1275–1288, 2012. [15] A. Cipitria, W. Wagermaier, P. Zaslansky, H. Schell, J. Reichert, P. Fratzl, D. Hutmacher, and G. Duda. BMPdelivery complements the guiding effect of scaffold architecture without altering bone microstructure in critical-sized long bone defects: a multiscale analysis.
Acta biomaterialia , 23:282–294, 2015.[16] P. G. Coelho, S. J. Hollister, C. L. Flanagan, and P. R. Fernandes. Bioresorbable scaffolds for bone tissue engineer-ing: optimal design, fabrication, mechanical testing and scale-size effects analysis.
Medical engineering & physics ,37(3):287–296, 2015.[17] M. R. Dias, J. M. Guedes, C. L. Flanagan, S. J. Hollister, and P. R. Fernandes. Optimization of scaffold design forbone tissue engineering: a computational and experimental study.
Medical engineering & physics , 36(4):448–457,2014.[18] J. Diestel and J. Uhl.
Vector Measures . American Mathematical Society, 1977.[19] P. Dondl, P. S. P. Poh, M. Rumpf, and S. Simon. Simultaneous elastic shape optimization for a domain splitting inbone tissue engineering.
Proc. A. , 475(2227):20180718, 17, 2019.[20] A. Ern and J.-L. Guermond.
Theory and practice of finite elements , volume 159. Springer Science & Business Media,2013.[21] P. Grisvard.
Elliptic problems in nonsmooth domains . SIAM, 2011.[22] K. Gr¨oger. A W ,p -estimate for solutions to mixed boundary value problems for second order elliptic differentialequations. Mathematische Annalen , 283(4):679–687, 1989.[23] J. K. Guest and J. H. Pr´evost. Optimizing multifunctional materials: design of microstructures for maximizedstiffness and fluid permeability.
International Journal of Solids and Structures , 43(22-23):7028–7047, 2006.[24] R. Haller-Dintelmann, C. Meyer, J. Rehberg, and A. Schiela. H¨older continuity and optimal control for nonsmoothelliptic problems.
Applied Mathematics and Optimization , 60(3):397–428, 2009.[25] H. Kang, C.-Y. Lin, and S. J. Hollister. Topology optimization of three dimensional tissue engineering scaffoldarchitectures for prescribed bulk modulus and diffusivity.
Structural and Multidisciplinary Optimization , 42(4):633–644, 2010.[26] M. Kassmann and W. Madych. Regularity for linear elliptic mixed boundary problems of second order. preprint ,2004.[27] D. H. Kempen, L. B. Creemers, J. Alblas, L. Lu, A. J. Verbout, M. J. Yaszemski, and W. J. Dhert. Growth factorinteractions in bone regeneration.
Tissue Engineering Part B: Reviews , 16(6):551–566, 2010.[28] V. Klika, M. A. P´erez, J. Garc´ıa-Aznar, F. Marˇs´ık, and M. Doblar´e. A coupled mechano-biochemical model for boneadaptation.
Journal of Mathematical Biology , 69(6):1383–1429, 2014.[29] O. A. Ladyzhenskaia, V. A. Solonnikov, and N. N. Ural’tseva.
Linear and quasi-linear equations of parabolic type ,volume 23. American Mathematical Soc., 1968.[30] C. Y. Lin, N. Kikuchi, and S. J. Hollister. A novel method for biomaterial scaffold internal architecture design tomatch bone elastic properties with desired porosity.
Journal of biomechanics , 37(5):623–636, 2004.[31] J. L. Lions.
Equations differentielles operationnelles: et probl´emes aux limites , volume 111. Springer-Verlag, 2013.[32] C. Marin, F. P. Luyten, B. Van der Schueren, G. Kerckhofs, and K. Vandamme. The impact of type 2 diabetes onbone fracture healing.
Frontiers in Endocrinology , 9:6, 2018.[33] L. A. Mills, S. A. Aitken, and A. H. R. Simpson. The risk of non-union per fracture: current myths and revisedfigures from a population of over 4 million adults.
Acta orthopaedica , 88(4):434–439, 2017.[34] A. Nauth, E. Schemitsch, B. Norris, Z. Nollin, and J. T. Watson. Critical-size bone defects: is there a consensus fordiagnosis and treatment?
Journal of orthopaedic trauma , 32:S7–S11, 2018.[35] M. Paris, A. G¨otz, I. Hettrich, C. M. Bidan, J. W. Dunlop, H. Razi, I. Zizak, D. W. Hutmacher, P. Fratzl, G. N.Duda, et al. Scaffold curvature-mediated novel biomineralization process originates a continuous soft tissue-to-boneinterface.
Acta biomaterialia , 60:64–80, 2017.[36] A. Petersen, A. Princ, G. Korus, A. Ellinghaus, H. Leemhuis, A. Herrera, A. Klaum¨unzer, S. Schreivogel,A. Woloszyk, K. Schmidt-Bleek, et al. A biomaterial with a channel-like pore architecture induces endochondralhealing of bone defects.
Nature communications , 9(1):1–16, 2018.[37] C. Pitt, F. Chasalow, Y. Hibionada, D. Klimas, and A. Schindler. Aliphatic polyesters. i. the degradation of poly( (cid:15) -caprolactone) in vivo.
Journal of applied polymer science , 26(11):3779–3787, 1981.[38] A.-M. Pobloth, S. Checa, H. Razi, A. Petersen, J. C. Weaver, K. Schmidt-Bleek, M. Windolf, A. ´A. Tatai, C. P.Roth, K.-D. Schaser, et al. Mechanobiologically optimized 3d titanium-mesh scaffolds enhance bone regeneration incritical segmental defects in sheep.
Science translational medicine , 10(423), 2018.[39] P. S. Poh, D. Valainis, K. Bhattacharya, M. van Griensven, and P. Dondl. Optimization of bone scaffold porositydistributions.
Scientific Reports , 9(1):9170, 2019.[40] Y. Qin.
Analytic inequalities and their applications in PDEs . Springer, 2017.[41] J. C. Reichert, M. E. Wullschleger, A. Cipitria, J. Lienau, T. K. Cheng, M. A. Sch¨utz, G. N. Duda, U. N¨oth, J. Eulert,and D. W. Hutmacher. Custom-made composite scaffolds for segmental defect repair in long bones.
InternationalOrthopaedics , 35(8):1229–1236, 2011.[42] E. Roddy, M. R. DeBaun, A. Daoud-Gray, Y. P. Yang, and M. J. Gardner. Treatment of critical-sized bone defects:clinical and tissue engineering perspectives.
European Journal of Orthopaedic Surgery & Traumatology , 28(3):351–362, 2018.
N EFFICIENT MODEL FOR SCAFFOLD-MEDIATED BONE REGENERATION 29 [43] J. A. Sanz-Herrera, J. M. Garcia-Aznar, and M. Doblare. A mathematical model for bone tissue regeneration insidea specific type of scaffold.
Biomechanics and Modeling in Mechanobiology , 7(5):355–366, 2008.[44] G. Savar´e. Regularity and perturbation results for mixed second order elliptic problems.
Communications in PartialDifferential Equations , 22(5-6):869–899, 1997.[45] H. K. Schwyzer, J. Cordey, S. Brun, P. Matter, and S. M. Perren. Bone loss after internal fixation using plates,determination in humans using computed tomography. In S. M. Perren and E. Schneider, editors,
Biomechanics:Current Interdisciplinary Research: Selected proceedings of the Fourth Meeting of the European Society of Biome-chanics in collaboration with the European Society of Biomaterials, September 24–26, 1984, Davos, Switzerland ,pages 191–195. Springer Netherlands, Dordrecht, 1985.[46] G. Stampacchia. Contributi alla regolarizzazione delle soluzioni dei problemi al contorno per equazioni del secondoordine ellittiche.
Annali della Scuola Normale Superiore di Pisa-Classe di Scienze , 12(3):223–245, 1958.[47] S. Stewart. Fracture non-union: A review of clinical challenges and future research needs.
Malaysian orthopaedicjournal , 13(2):1, 2019.[48] T. Terjesen, A. Nordby, and V. Arnulf. Bone atrophy after plate fixation: Computed tomography of femoral shaftfractures.
Acta Orthopaedica Scandinavica , 56(5):416–418, 2009.[49] V. Viateau, G. Guillemin, V. Bousson, K. Oudina, D. Hannouche, L. Sedel, D. Logeart-Avramoglou, and H. Petite.Long-bone critical-size defects treated with tissue-engineered grafts: A study on sheep.
Journal of OrthopaedicResearch , 25(6):741–749, 2007.[50] X. Wang, S. Xu, S. Zhou, W. Xu, M. Leary, P. Choong, M. Qian, M. Brandt, and Y. M. Xie. Topological designand additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: A review.
Biomaterials ,83(c):127–141, Mar. 2016.[51] S. R. Yu, M. Burkhardt, M. Nowak, J. Ries, Z. Petr´aˇsek, S. Scholpp, P. Schwille, and M. Brand. FGF8 morphogengradient forms by a source-sink mechanism with freely diffusing molecules.
Nature , 461(7263):533–536, 2009.[52] G. Zimmermann and A. Moghaddam. Trauma: non-union: new trends. In
European instructional lectures , pages15–19. Springer, 2010.(Patrick Dondl)
Abteilung f¨ur Angewandte Mathematik, Albert-Ludwigs-Universit¨at Freiburg, Hermann-Herder-Strasse 10, 79104 Freiburg i. Br.
Email address : [email protected] URL : https://aam.uni-freiburg.de/agdo/ (Patrina S. P. Poh) Julius Wolff Institute for Biomechanics and Musculoskeletal Regeneration, Charit´e –Univerist¨atsmedizin Berlin, Berlin, Germany
Email address : [email protected] URL : https://jwi.charite.de/metas/person/person/address detail/poh/ (Marius Zeinhofer) Abteilung f¨ur Angewandte Mathematik, Albert-Ludwigs-Universit¨at Freiburg, Hermann-Herder-Strasse 10, 79104 Freiburg i. Br.
Email address ::