Governing Equations of Tissue Modelling and Remodelling: A Unified Generalised Description of Surface and Bulk Balance

Several biological tissues undergo changes in their geometry and in their bulk material properties by modelling and remodelling processes. Modelling synthesises tissue in some regions and removes tissue in others. Remodelling overwrites old tissue material properties with newly formed, immature tissue properties. As a result, tissues are made up of different “patches”, i.e., adjacent tissue regions of different ages and different material properties, within evolving boundaries. In this paper, generalised equations governing the spatio-temporal evolution of such tissues are developed within the continuum model. These equations take into account nonconservative, discontinuous surface mass balance due to creation and destruction of material at moving interfaces, and bulk balance due to tissue maturation. These equations make it possible to model patchy tissue states and their evolution without explicitly maintaining a record of when/where resorption and formation processes occurred. The time evolution of spatially averaged tissue properties is derived systematically by integration. These spatially-averaged equations cannot be written in closed form as they retain traces that tissue destruction is localised at tissue boundaries. The formalism developed in this paper is applied to bone tissues, which exhibit strong material heterogeneities due to their slow mineralisation and remodelling processes. Evolution equations are proposed in particular for osteocyte density and bone mineral density. Effective average equations for bone mineral density (BMD) and tissue mineral density (TMD) are derived using a mean-field approximation. The error made by this approximation when remodelling patchy tissue is investigated. The specific signatures of the time evolution of BMD or TMD during remodelling events are exhibited. These signatures may provide a way to detect remodelling events at lower, unseen spatial resolutions from microCT scans.


Introduction
Tissue growth, renewal, and shape adaptation are common traits to many biological tissues and biomaterials. These traits are enabled by the processes of tissue modelling (tissue generation or destruction) and tissue remodelling (renewal by coordinated destruction and regeneration).
Tissue growth enables us to be born small and to grow to maturity [1]. Tissue shape adaptation and renewal enables structural reorganisation, maturation, and self-repair, which are important factors of tissue function. For example, bone tissues adapt their shape and microstructure to mechanical loads to offer strength with minimal weight, and they repair microcracks to prevent structural damage. Muscles and tendons adapt their mass and fibre structure to the forces they transmit [2,3]. Extracellular matrix (ECM) modelling and remodelling helps cells to migrate [4] and it give cells control over local stress fields, for example to provide stress shielding [5]. Modelling and remodelling are often associated with the evolution of internal or external tissue boundaries (Fig 1), such as in wound repair and reconstruction of damaged ECM, which proceed as self-organised wave propagations [6,7]. Cancer invasion breaks down normal tissues boundaries, rearranging their architecture and affecting their function.
While some tissues are renewed in a linear fashion with creation consistently occurring in one region and removal occurring in another (e.g., nail, hair, skin), other tissues have more complex patterns of creation and removal (e.g. ECM, bone), resulting in tissue heterogeneities that reflect the history of their generation.
The evolution of tissue material properties is challenging to grasp within a single mathematical modelling framework due to tissue heterogeneities and moving boundaries. The record of maturing tissue properties may suddenly and locally be erased and overwritten with immature tissue material, creating internal discontinuities in bulk material properties within the tissue. Ordinary differential equations (ODEs) describe the time evolution of spatially averaged tissue properties, but it is unclear how changes occurring at boundaries are reflected in such spatial averages. Partial differential equations (PDEs) describe the spatio-temporal evolution of tissue properties. However, to represent discontinuities at moving interfaces, these equations must possess singular terms. The nature of these singularities is the main topic of this paper. Mathematical and computational models typically avoid such singularities by resorting to (i) volume of fluid methods or mixture theory, which represent the evolution of continuous partial fractions that in effect smooth out boundaries; or (ii) discrete models, for which discontinuities pose no particular problem [8][9][10][11].
In this paper, general governing equations are proposed to describe the evolution of tissue geometry and tissue properties through bulk maturation processes and through formation and resorption processes localised at tissue boundaries. The novelty of these equations is in accounting for nonconservative, discontinuous surface balance due to creation and destruction of quantities at moving boundaries. The discontinuities associated with tissue modifications at boundaries are captured by singular terms, namely, surface distributions, to be understood in the sense of distribution theory [12][13][14]. These generalised material balance equations are widely applicable and as many biochemical and transport processes as necessary can be included for a particular application. The formalism developed in this paper is anticipated to find particularly useful applications in tissue engineering, biofabrication, and investigations of bioscaffold integration and remodelling [15][16][17].
Surface distributions have been introduced in the context of non-equilibrium thermodynamics, interfacial conservation equations, and Stefan problems with several different representations and degrees of rigour, which has sometimes led to confusion between authors [13,[18][19][20][21][22][23][24]. Important properties of the surface distribution and the equivalence of these representations are shown in Appendix A.
A main advantage of formulating governing equations of tissue modelling and remodelling over discrete models, is that these equations lend themselves to mathematical analysis. We will see that this formalism enables the systematic derivation of equations governing the temporal evolution of spatial averages of tissue properties, such as density in a representative elementary volume of tissue, as well as total tissue volume. This derivation reveals that traces that tissue removal is localised at the tissue boundary are retained in the resulting ODEs, preventing these equations from being written in closed form. The error made by closing these equations with a mean-field approximation is investigated.
A concrete application of this formalism is developed to describe the evolution of bone tissues under bone modelling and bone remodelling processes, with a focus on two applications of particular interest in bone: (i) evolution equations of bone-embedded cells (osteocytes); and (ii) bone mineralisation. The example (i) extends the model of osteocyte formation and viability introduced in Ref. [14] by including the effect of local removal by bone resorption. This extension enables the representation of heterogenous bone states and their evolution during bone remodelling. The example (ii) is particularly important as experimental and clinical bone scans typically provide a measure of bone mineral density, averaged over spatial regions corresponding to the scanners' resolution [25].

Material balance of local tissue properties
Consider a local material property of the tissue η(r, t) at position r in space and at time t. The value of η is assumed to be zero out of the tissue's spatial extent. Conceptually, the definition of η(r, t) involves small representative elementary volumes within which the material property is spatially averaged. These volumes are assumed large enough to contain many molecules so that the property is well-defined, but small enough so that spatial inhomogeneities occurring at a larger scale when r is varied are not averaged out. The continuum model formally takes the limit to zero of these elementary volumes to define η(r, t) at every point r of the continuous space [26][27][28]. In this limit, the molecular detail is omitted and properties such as η(r, t) become generalised functions governed by equations to be interpreted in the sense of distributions [13,29,30].
Tissue modelling, remodelling, and maturation modify η(r, t) in several ways. Away from tissue boundaries, changes in η are due to bulk processes such as chemical reactions and internal transport phenomena. Most of a tissue's heterogeneity is not due to such processes. It is due to the dynamic nature of tissue synthesis. Different regions of the tissue are created at different times. They have different properties η reflecting different ages and different biological contexts at creation. Tissue heterogeneities seen in a property η are a record of when, where, and how the tissue was synthesised or modified. We will assume here that this synthesis or modification occurs by cellular action located at a boundary S(t), which may be an internal boundary within the tissue, or the tissue boundary (Fig 1). We assume quite generally that cellular action at S(t) sets a new value of η there (see Fig 2). The equation that governs this process is given by:

@ @t
Zðr; tÞ ¼ DZðr; tÞ vðr; tÞd SðtÞ ðrÞ; ð1Þ where Δη is the change in η occurring at S(t) by the cellular action, v is the normal velocity of S (t), and δ S(t) is the surface distribution, formally zero everywhere except at S(t), where it is infinite. This singularity indicates the discontinuous nature of η at S(t). It will be responsible for the creation of sharp internal boundaries within the tissue when the normal velocity of S(t) changes sign, for example at reversals between tissue resorption and tissue formation. Mathematically, δ S (r) is a distribution defined such that it maps any test function φ(r) to the real value given by the surface integral of φ over S [13,[18][19][20]: Several properties of δ S are presented in Appendix A. The justification of Eq (1) is given by the following jump property: Jump property. Let η be governed by Eq (1) and let t Ã = t Ã (r) be the arrival time of the boundary S(t) at point r, such that r 2 S(t Ã ) (see Fig 2). Then, the value of η at r is constant except at time t Ã where it jumps by the quantity Δη: where t Ã AE ¼ t Ã AE , and ! > 0. The tissue material property η(r, t) at r jumps by the quantity Δη at the arrival time t* (r). The sign of the jump depends on the direction of propagation of the boundary. Right: in a small neighbourhood of r, it is possible to choose a local coordinate system with components parallel to and perpendicular to S(t*(r)) such that the infinitesimal volume element dr = dσ|v|dt (see also explanations in the text).
This property is demonstrated in Appendix A 'Jump property'. One can see from Fig 2 that the sign of the jump in η depends on the direction of propagation of the boundary, and therefore on the sign of v: if the boundary in Fig 2 travels toward the right, the value of η at r increases by Δη > 0 at the passage of the boundary; if the boundary travels toward the left, the value of η at r decreases by Δη at the passage of the boundary. In practice, the normal velocity of cells at tissue boundaries can always be assumed positive and the sign of the jump is then solely determined by that of Δη: a positive sign represents formation, a negative sign represents resorption. The jump in η at r only depends on the values of v and Δη prevailing at the arrival time t Ã . At any other time than the arrival time t Ã , Eq (1) implies that η at r is constant.
The full balance of a general tissue property η is obtained by adding to Eq (1) further conservative and nonconservative processes that modify η at other times. For illustration, consider a tissue occupying a region O(t) in space with boundary S(t) = @O(t). The tissue is assumed to change shape due to tissue formation and resorption occurring at specific regions of S(t). It is also assumed to change its material properties due to maturation. The surface and bulk balance of a property η of this tissue can be formulated based on Eq (1) as follows: • Tissue formation. New tissue is deposited on S(t) with a normal velocity v = v f > 0 and an initial material property η f (r, t). Both v f and η f are determined by the synthesis process (e.g., cell secretion). By Eq (1), the rate of change in η due to this process is η f v f δ S(t) .
• Tissue resorption. Existing tissue is resorbed from S(t) with a normal velocity v = −v r < 0 determined by the removal process (e.g., cell-driven chemical dissolution or mechanical wear). The property η drops from its current value to zero. By Eq (1), the rate of change in η due to this process is η due to this process is −η − v r δ S(t) , where η − is the value of η probed at an infinitesimal inward offset of S(t).
• Tissue maturation. After new tissue synthesis, η evolves according to biochemical and mechanical processes specific to η, until it is removed by resorption. The rate of change in η due to this process is assumed to be given by a maturation law F ðZ; r; tÞ.
The evolution of η is given by summing up these contributions: Some regions of S(t) may undergo formation while others may undergo resorption simultaneously. Since these regions may not overlap, the normal velocity is given everywhere by where v f and v r correspond to the positive and negative parts of v. The evolution of the tissue's shape is univocally determined by v [31,32]. The regularisation η − in the resorption term is necessary because η is discontinuous at S(t). It ensures that the value of η to remove during resorption is probed at a point lying just within the tissue rather than where it jumps to 0. In the sequel, we will omit this regularisation from the notation with the convention that η takes the value η − whenever it is evaluated at S(t).
The maturation law F represents a general bulk balance of η, which may include nonconservative processes such as chemical reactions, and conservative processes due to transport phenomena within the tissue.
Nonconservative vs conservative surface balance. Eq (4) is a generalised balance equation that explicitly accounts for nonconservative processes occuring at moving interfaces due to creation and destruction of material. It has similar surface terms as conservation equations in multiphase systems and Stefan problems [18-20, 24, 33-37]. The main difference is that surface terms in these systems are inherently conservative. They represent jump conditions necessary to enforce mass conservation at the interface. To illustrate the difference, consider the general transport theorem that expresses the total variation of η in an evolving domain O(t) [35]. Taking O(t) to follow the material velocity of η such that there is no influx or outflux of η through @O(t), one has d dt where the term in the right hand side represents change in η within O(t) due to nonconservative phenomena such as chemical reactions. The surface integral can be rewritten is an arbitrary region of the substance, η must be governed locally by: In Eq (6), the conservative balance of η imposes the fact that the jump in η at locations of S(t) with a positive normal velocity v = v f , is the value of η at an infinitesimal inward offset of @O (t), rather than an independent value η f determined by nonconservative processes as in Eq (4). Furthermore, the normal velocity v of the boundary in Eqs (5)-(6) is determined by the material velocity of the substance η, whereas in Eq (4), it is determined by the independent processes of new tissue formation and resorption occuring at the interface. Naturally, both conservative and nonconservative surface balance terms may in general be present in the balance of a property.

Evolution of spatially averaged tissue properties
Many mathematical models describe the evolution of tissue properties in time only. These models implicitly assume that the property is distributed homogeneously in the tissue. Eq (4) enables us to derive systematically the time evolution of spatial averages of patchy tissue properties, and to investigate the error made by assuming tissue homogeneity. (See Refs [21,22,33,35,38] for volume averaging theorems in conservation equations.) Let V be a fixed mesoscopic or macroscopic representative elementary volume and O(t) be the volume occupied by the tissue in V, with boundary S(t) = @O(t). The tissue volume fraction in V is (We use V, O(t), and S(t) to denote both the region in space and the measures |V|, |O(t)|, and | S(t)| for simplicity.) Two spatial averages of η can be defined based on V and O(t): The average hηi V may integrate η over regions devoid of tissue, where η = 0. It is thus related to hηi O through the tissue volume fraction: Differentiating Eq (8) with respect to t and using Eq (4) gives: where hÁi S ¼ 1 SðtÞ R SðtÞ ds Á is the average value over S(t). The surface density S(t)/V (also called specific surface) is an important characteristic of porous media. For example in bone tissues, it is related to the propensity to remodel [39,40]. Let S f (t) and S r (t) denote the forming and resorbing surfaces of S(t), i.e., the portions of S(t) at which v f 6 ¼ 0 and v r 6 ¼ 0, respectively. Eq (11) can be rewritten as: If η is taken to be the indicator function 1 OðtÞ of O(t), then h1 OðtÞ i V ¼ f ðtÞ and Eqs (11), (12), together with the balance equation of the indicator function (Eq (50) in Appendix A 'Balance equation of the indicator function of an evolving domain') determine the evolution of the tissue volume fraction f(t): Eq (13) shows in particular that the tissue volume O(t) = Vf(t) evolves according to: i.e., tissue volume changes at a rate equal to the tissue surface area multiplied by the average normal velocity, as expected. To determine the evolution of averages defined with O(t) as volume referent, note that from Eq (10): Using Eqs (12) and (13), one obtains Eqs (12) and (16) show that the evolution of spatial averages of patchy tissues cannot be written in closed form even when F ðZÞ is linear, i.e., even when hF ðZÞi ¼ F ðhZiÞ. Indeed, due to resorption, changes in hηi V or hηi O depend on the value of η deposited last, occurring in the factor hv r Zi S r , rather than on the current volume average. This hysteresis of the evolution of averages is due to the fact that tissue resorption proceeds from the tissue surface, and thus removes a value of η that depends on when and how it was first deposited. In Section 'Meanfield approximation of bone mineral density', the error committed when closing the equations using a mean-field approximation is studied on bone mineral density. Note that if the material property does not mature (F 0) and if a constant value η f is generated during tissue formation, then there is no hysteresis, and trivially,

Application to bone tissue
Bone is a dynamic tissue that sustains lifelong changes in its microstructure and in its material properties [41]. At the cellular scale, bone is composed of (i) bone matrix, infiltrated with minerals and with the osteocyte network; and (ii) vascular pores, containing soft tissues and cells. Changes in bone microstructure occur by dissolution of old bone matrix by bone-resorbing cells (osteoclasts) and deposition of new bone matrix by bone-forming cells (osteoblasts) [4,41,42]. Changes in material properties of newly deposited bone occur by matrix maturation such as collagen fiber re-arrangement, mineralisation, accumulation of micro-cracks, and maturation of osteocytes [41,42].
Bone remodelling turns over bone tissue slowly, at rates of 5-30%/year. This allows bone matrix to undergo significant changes in material properties before being renewed. As a result, the state of bone is "patchy": it contains many internal boundaries separating tissue regions of different ages, which reflect the history of their formation and resorption processes. These different tissue regions are called bone structural units or osteons [41][42][43].
During bone modelling and remodelling, the bone surface S(t) between bone matrix and vascular pores evolves by the action of osteoblasts and osteoclasts. The normal velocity of S(t) is given by where ρ Ob (r, t) and ρ Oc (r, t) are the surface density of osteoblasts and osteoclasts (number per unit surface), k f (r, t) is the secretory rate (volume formed per osteoblast per unit time), and k r (r, t) is the resorption rate (volume dissolved per osteoclast per unit time) [14].
Osteocyte density. Osteocytes are tissue-embedded cells believed to sense and transduce mechanical strains of bone matrix to osteoblasts and osteoclasts. Osteocytes reside in small cavities and channels within bone matrix, making up a porosity of about 1-2% [44]. No modelling or remodelling is initiated at these micropore surfaces. Osteocytes are generated with new bone matrix during formation. They can be viewed as a bone material property, generated initially with density Ot f . The spatio-temporal evolution of osteocyte density Ot (N.Ot/BV in bone histomorphometric standards [45]) is governed by: The last term accounts for apoptosis (cell death) occurring with rate A(r, t). In Ref. [14], a similar evolution equation for osteocyte density was proposed, but no resorption was accounted for. The first term was modelled as D burial ρ Ob δ S(t) to represent the fact that osteocytes are osteoblasts that become buried during bone formation, where D burial (r, t) is the burial rate, i.e., the probability per unit time for an osteoblast to become trapped in bone as an osteocyte. By identification with the term Ot f v f δ S(t) in Eq (18), one immediately finds that the density of osteocyte generated at the moving deposition front is given by as obtained in [14]. Eq (19) holds generally for the density of any inclusion deposited by osteoblasts in bone matrix at rate D burial , and by extension, for any inclusion in tissue or material synthesised at an interface. This density does not explicitly depend on surface curvature and osteoblast density. This is particularly relevant for biological tissues and biomaterials, since consistent inclusion densities can be generated in complex geometries and nonconstant populations of tissue-synthesising cells simply by maintaining the cell-specific properties D burial and k f constant.
Eq (18) was solved numerically in one spatial dimension (z) when A and Ot f are constant, and v = v(t) oscillates between two values (Fig 3). In where t Ã (r) is the arrival time at r, BV(t) is the spatial region occupied by bone at time t, and 1 BVðtÞ is the indicator function of BV(t). This solution was derived in Ref. [14] like Eqs (21)- (25) for bone mineral density below. In Fig 3b, v(t) oscillates between a positive and a negative value: there is an alternation of bone tissue formation and bone tissue resorption. Resorption introduces sharp discontinuities in osteocyte density in adjacent regions, resulting in a bone matrix composed of distinct tissue layers ("patches"). These patches are due to the fact that tissue lying under resorbing surfaces keep maturing. When resorption stops and new tissue forms, there is an age gap between the underlying tissue and new tissue.
The analytic solution Eq (20) holds within each patch region, which may shrink during resorption. The solution in the whole space can be constructed piecewise. However, this requires book-keeping of the time and locations at which there is reversal between resorption and formation to identify patches. Such book-keeping is tedious and impractical in higher dimensions as tissue formation events may be generated at different times and locations of the surface. The governing Eq (18) can represent these patches without explicitly needing the information of resorption-formation reversals. It can also handle more elaborate situations, such as nonlinearities and complex couplings. Bone mineral density. New bone is formed initially as an unmineralised collagen matrix. This unmineralised matrix matures and gradually incorporates minerals to become hard bone tissue. Mineralisation first increases rapidly due to the deposition of mineral pellets by cells during formation. It then continues to increase over much larger time scales by crystal growth [42]. Mineral density is an important bone material property. It is measured clinically as an indicator of skeletal integrity, for example in osteoporosis [25]. Assuming that new bone tissue is infiltrated with an initial density of mineral pellets m f (r, t), the spatio-temporal evolution of bone mineral density is governed by: The mineralisation law F miner determines the evolution of mineral density after the initial pellet deposition. Without resorption, m(r, t) is solution of the initial value problem @ @t mðr; tÞ ¼ F miner ðmÞ; 8t > t Ã ðrÞ; ð22Þ The initial value Eq (23) expresses the jump property Eq (3) at time t = t Ã (r) due to the surface balance term m f v f δ S(t) : at t = t Ã (r), m jumps from 0 to m f . After the initial mineral deposition, we assume that bone mineral density increase until it reaches a maximum mineral density m max . We model this mineralisation process by exponential saturation: If m max (r) is independent of time in Eq (24), the solution to Eqs (22)-(24) is: where m f is evaluated at m f (r, t Ã (r)). In reality, m max is likely to be a function of time. It is believed to be regulated by osteocytes and their dendritic processes [42,46]. Experimental determinations of the increase in mineral density with time in newly deposited bone tissue exhibit two time scales (Fig 4) [47]. While explicit fitting functions for t 7 ! m (r, t) have been proposed with great accuracy to experimental data [47], these fitting functions do not satisfy a simple mineralisation kinetics law of the type Eq (22). We assume instead that mineralisation is described by the exponential saturation law Eq (24) with distinct characteristic times at these two time scales. The constants m f , m max , and k m in Eq (24)   . Bone mineral density was evolved using the short-time mineralisation parameters. Bone-resorbing cells first create a cavity in a mineralising bone tissue substrate. Bone-forming cells then deposit new tissue. The new tissue contrasts with the older substrate by its lower mineral content. After remodelling has completed, the tissue is clearly made up of two distinct patches. Within each patch, the mineral density keeps increasing and is continuous, but it is discontinuous at the line corresponding to the deepest location reached by resorption. In bone, this line of reversal between resorption and formation is called the cement line. The patch of newly formed bone is called a secondary osteon, or bone structural unit [42].
The overall bone balance after the remodelling event in Fig 5 is approximately zero. However, the interface has changed. Small changes in the interface are likely to occur in bone remodelling even without bone loss or gain. Indeed, bone remodelling is regulated by several processes of biochemical, geometrical, and mechanical nature, which affect the generation and   [47] with contour lines shown every 1 Ca wt%. The color scale of the interface is the normal velocity, normalised by the maximum absolute value in this simulation. At t = 50 days, remodelling is initiated with osteoclasts (red) starting to resorb bone matrix until t = 90 days. At t = 70 days, osteoblasts (blue) are activated towards the rear and start refilling the resorbed cavity until t = 120 days, at which point the interface is still and remodelling has completed. Because bone tissue resorption has removed a portion of bone adjacent to mineralising tissue, newly formed bone contrasts by its mineral content with surrounding tissue. The end state of the bone matrix is made up of two distinct patches, within which mineral density is continuous. coupling of bone-resorbing and bone-forming cells. These regulatory processes were not modelled here, see Refs [48][49][50][51] for more biologically accurate mathematical models of cell population dynamics in BMUs. Fig 6 shows a portion of bone that underwent two bone remodelling events in twenty years, roughly corresponding to a turnover rate of 10%/year [41,42]. The cell populations were assigned so as to emulate remodelling events without net bone gain or loss (Appendix B). Bone mineral density was evolved using the long-time mineralisation parameters. After the second remodelling event, the tissue is made up of three distinct patches: the old bone substrate, bone renewed by the first remodelling event, and bone renewed by the second remodelling event. Part of the bone renewed by the first remodelling event was removed and replaced by newer bone during the second remodelling event. This kind of variegated state of bone matrix is typical, as observed by microradiographs [52][53][54], quantitative back-scattering electron microscopy [47,55], and micro-computed tomography [56,57]. Quantities recorded in bone during formation are gradually overwritten with newer content. This constitutes a loss of information: osteocyte density for example records the ratio of burial rate to secretory rate that is current at the time of formation, see Eq (19) [14]. On the other hand, bone tissue patches provide other information such as the age or turnover rate of the tissue.
As in the one-dimensional simulation, the governing Eq (21) can represent tissue patches without needing the information of the time and locations of resorption-formation reversals. Book-keeping patch location in space and time is particularly complicated in situations such as

Bone tissue spatial averages
Eqs (12)- (16) are valid in general. Here we specialise them to bone tissue using Eq (17) and further assume that the secretory rate k f and dissolution rate k r are constant. We also denote V by TV (tissue volume) and O(t) by BV (bone volume) to follow bone histomorphometric conventions [45]. Under these assumptions: where N.Ob and N.Oc are the number of osteoblasts and osteoclasts in TV. Using Eq (26) in Eq (13), bone volume fraction evolves as: Eq (27) provides a microscopic justification of the equation df ðtÞ dt ¼ k f Ob À k r Oc used in the literature, where Ob and Oc are average cell densities in a representative elementary volume [58][59][60][61].
Mean-field approximation of bone mineral density. Current conventional microCT scanners have millimetric to submillimetric resolution. They effectively measure local spatial averages of bone mineral densities. If soft tissues are included in the average, measurements refer to 'bone mineral density' (BMD). If soft tissues are excluded, measurements refer to 'tissue mineral density' (TMD) [62]. Thus: The evolution of hmi TV and hmi BV by Eqs (12) and (16) depends on the patchy state of bone, and so on remodelling history. However, if bone mineral density is not too inhomogeneous in BV, we can close Eqs (12) and (16) by making the mean-field approximation: With the mineralisation model used in Figs 4-6, which assumes k f , k r , k m , and m f constant, we have: where the last equalities in each line used Eq (26) and the first equality in the second line used the mean-field approximation Eq (29). Eq (12) thus becomes: and Eq (16) becomes: For given average densities of osteoblasts and osteoclasts, f is given by Eq (27), and Eqs (30), (31) are now self-consistent. The time evolution of hmi BV and hmi TV found by explicitly averaging the microscopic models (21), (24) exhibits specific model elements in the different terms of Eqs (30) and (31). These elements could easily be missed when heuristically formulating a temporal model directly: 1. The factor 1/f in Eq (30) is due to Eq (10); 2. The factor f multiplying m max in Eq (30) is due to the fact that F miner ðmÞ is not linear in m; it is dicontinuous at m = 0. In fact, F miner is such that m max = 0 out of BV, so that hm max i TV = fm max . (31) is independent of resorption. The dependence on formation corresponds to the relaxation of hmi BV towards the value deposited m f . The relaxation rate is proportional to the bone formation rate k f hObi BV and to 1/f. The lower the bone volume fraction f, the quicker it is to replace the current average hmi BV with new values m f . Except during the remodelling events, hmi TV and hmi BV increase due to the mineralisation law Eq (24). The large dips in hmi TV are due to the changes in f during the remodelling events. Both hmi TV and hmi BV have decreased values at the end of each remodelling event compared to the value prior to the remodelling event. This is due to the presence of new, lower-mineralised bone after remodelling. The numerical solutions of the mean-field Eqs (30) and (31) (scaled by either f or 1/f according to Eq (10)) are undistinguishable. The mean-field solutions differ from the averaged spatio-temporal solution by 0.01% at the first remodelling event, and by 0.04% at the start of the second remodelling event (see insets). These differences are attributed to the different numerical integrations required by the solutions. However, if model elements listed in points (i)-(iii) above are missed, the mean-field solutions can differ dramatically (not shown).

The evolution of hmi BV in Eq
Just before the first remodelling event, bone mineral density is approximately homogenous. When this bone is remodelled, there is no qualitative difference between the averaged spatiotemporal solution and the mean-field approximations (top-left inset). However, just before the second remodelling event, bone mineral density is distributed heterogenously across two patches (Fig 6a). When this bone is remodelled, the average mineral density at the surface is significantly lower than the bone volume average. The removal of bone near the surface during resorption thus accelerates the increase in hmi BV in a first stage (Fig 7b, solid line in bottomright inset), before hmi BV decreases due to new, lower-mineralised bone being deposited during formation. This initial accelerated increase in hmi BV is missed by the mean-field approximations (interrupted lines).
The behaviours of BMD and TMD around t = 10yr and t = 20yr in Fig 7 represent typical signatures of isolated remodelling events that could in principle be seen by in-vivo microCT of resolution TV. These time signatures may therefore provide a way to detect remodelling events that occur at lower, microscopic spatial resolutions not seen in the scans. Current in-vivo technologies remain limited in the number of timepoints and the transient behaviours during the remodelling events may be missed. However, sawtooth-like changes in mineral density may still be detected. This would require a resolution-accurate co-registration of scans taken at different timepoints. Note that concurring remodelling events within a voxel TV could smear out individual remodelling signatures.

Conclusions
This paper shows that the evolution of tissue geometry and tissue material properties under modelling and remodelling processes can be captured by a single, general mathematical framework. Tissue heterogeneities due to different tissue ages and different biological contexts at creation are represented in this framework by functions of space and time with discontinuities at internal or external boundaries. The equations governing the evolution of these functions are singular differential equations, in the sense of distribution theory. The surface distribution δ S and associated jump property Eq (3) enable a modular approach to formulating governing equations of complex tissues and biomaterials. They extend conventional balance equations with nonconservative processes localised at moving boundaries. This enables 'continuum model' notations to be employed despite the occurrence of discontinuities at surfaces, much like the Dirac distribution enables continuum notations to be employed in discrete systems [29].
Internal tissue boundaries that separate regions generated at different times are created at reversals between resorption and formation. These boundaries arise naturally from the governing Eq (4). In contrast, analytical solutions require to book-keep the time and location of these reversals to construct the solution piecewise from continuous patches.
A distinction is sometimes made in biology between tissue modelling and tissue remodelling. From the point of view of following the evolution of tissue properties, these processes do not need to be distinguished, so long as the effect of removal and formation on tissue properties are identical in both situations. This is the case of bone tissues, for which remodelling can be seen as a coordinated sequence of small resorptive and formative modelling processes [41]. The evolution of bone tissue during modelling or remodelling is thus described mathematically by the same set of equations, the difference being in the timing and location of the resorption and formation processes. Here, these were assumed given. In practice, this information may come from experimental data, or from further mathematical models of the populations of bone-resorbing and bone-forming cells. Generally, the governing Eq (4) needs to be supplemented with information on the specific processes involved in the formation and resorption kinetics of the tissue, which determine the normal velocities v f and v r , and the value η f of newly formed tissues.
Exact governing equations for the evolution of spatial averages of the tissue were obtained by integrating the spatio-temporal Eq (4). These average equations are not self-consistent due to the heterogeneous nature of the tissue, but can be closed by a mean-field approximation such as Eq (29). The degree to which the mean-field approximation is well satisfied depends on the degree of inhomogeneity of the tissue. Caution should be exercised whenever a tissue property changes over time scales that are faster than typical remodelling rates, which results in patchy states such as in Figs 5, 6. While the discrepancy due to the mean-field approximation in Fig 7b is small, such discrepancies would accumulate with further remodelling events. Fig 7 is a prediction of the type of BMD or TMD signatures that could be detected by in-vivo microCT scans when the bone undergoes remodelling at lower, unseen length scales.
The spatial and temporal scales at which the formalism presented in this paper is valid depend on the adequacy of the continuum model to represent a particular application at these scales. The bone tissue examples presented here were considering boundaries to be the bone-vascular interface. At a lower scale, boundaries may represent the secretory areas of a cell's membrane. At a higher scale, boundaries may represent the overall shape of an organ. This formalism is applicable to many other systems in which a material is created and destroyed from its surfaces while undergoing changes in the bulk. This includes tissues and biomaterials such as ECM remodelling, tooth development, the generation and biomineralisation of shells, bioscaffolds, but also non-biological systems, such as sedimentation, 3D printing, etching, and chemical adsorption.

Appendix A: Properties of the surface distribution
This appendix presents a few properties of the surface distribution δ S defined by Eq (2). See Refs [18][19][20] and Sec. 8.4 in Ref [13]. Intuitively, the surface distribution is similar to the Dirac distribution except that it is formally infinite on a N − 1 manifold embedded in R N . Such a manifold is usually called a 'hypersurface'. We will refer to it as a 'surface' for simplicity. It corresponds to a curve when N = 2 and a point when N = 1. Integrating the surface distribution δ S over N-dimensional space with a test function only retains the function's values on the N − 1 dimensional surface, and integrates these values with respect to the measure defining N − 1 dimensional area [63]. It is important to contrast the surface distribution with the Dirac distribution, which in all dimensions returns the value of a test function at a single point. We refrain from using the terminology 'Dirac' to refer to the surface distribution δ S to avoid potential confusion. In Section 'Local curvilinear partition of space', we first demonstrate the local curvilinear partition of space, Eq (35). This relation means that in effect, integrating over space with δ S removes spatial components normal to the surface S, see Eq (43). In Section 'Representations of the surface distribution', we mention several different representations of the surface distribution found in the literature.

Jump property
We first demonstrate the jump property enunciated in Eq (3).
Proof. In one spatial dimension, the interface is a point of coordinate S(t). Eq (3) is obtained by integrating Eq (1) over t 2 ½t Ã À ; t Ã þ and by using d SðtÞ ðxÞ ¼ dðx À SðtÞÞ ¼ 1 jvj dðt À t Ã Þ, where |v| = |S 0 (t Ã )|. To prove the jump property in higher dimensions, we first replace the running time variable t in the right hand side of Eq (1) by t Ã : only the values of Δη and v at the arrival time t Ã contribute to the change in η at r. Indeed, for any function φ(r, t): Z drd SðtÞ ðrÞ φðr; tÞ ¼ where the second equality in Eq (32) uses the fact that any point r 2 S(t) has the arrival time t Ã (r) = t. The pointwise notation in Eq (1) is elucidated in the theory of distributions by integrating over space with a smooth kernel function δ n of unit integral, and of support tending to the single point {r} as n ! 1. The sequence δ n is called a regular sequence converging to Dirac's delta distribution [13]: lim n ! 1δ n (r 0 − r) = (r 0 − r). The meaning of Eq (1) is thus @ @t Zðr; tÞ ¼ lim To calculate the jump in η induced by the passage of S(t) through r at t = t Ã , we integrate Eq (33) over t 2 ½t Ã À ; t Ã þ and use the definition Eq (2) The final step consists in partitioning space in a neighbourhood of r by a set of parallel and perpendicular coordinates to the interface S(t Ã ), such that This partioning is visually intuitive (see Fig 2). It is proved in Sec. 'Local curvilinear partition of space'. With Eq (35), we finally obtain where V (r, t Ã ) corresponds to the region in space swept by S(t) during t 2 ½t Ã À ; t Ã þ . If t Ã is not the arrival time at r, then for sufficiently small and sufficiently large n, the support of δ n(r 0 − r) is not contained in V (r, t Ã ), the integral in the right hand side of Eq (36) is zero, and η is unchanged. An alternative derivation of the jump property Eq (3)

Local curvilinear partition of space
We first show that it is possible to define a local curvilinear coordinate system around the point r with N − 1 coordinates parallel to S(t) and one coordinate perpendicular to S(t) for t around the arrival time t Ã at r. This local curvilinear coordinate system defines a local partition of the space around r such that an infinitesimal volume element dr will be represented by |v|dtdσ.
Let ψ(u, t) be a local parameterisation of the manifold S(t) around the point r, where u belongs to an open subset V & R NÀ1 . Under appropriate regularity conditions on the normal velocity v and on S(t) it is always possible to choose the time dependence of ψ such that the curves t 7 ! ψ(u, t) define trajectories normal to S(t) around t Ã for all u, by solving the differential equation @c @t ¼ vðc; tÞnðc; tÞ ð 37Þ from an initial parameterisation. In particular, S(t) must have no 'corners' in a small neighbourhood of r, S(t) must be an 'evolving hypersurface' around r [64]. The parameterisation thus obtained, can be seen as a coordinate transformation that maps the curvilinear coordinates (u, t) to the cartesian coordinates r 0 . In doing so, time lines become replaced by the distance travelled along trajectories perpendicular to S(t) and lines parameterised by u i become replaced by the distance travelled along trajectories parallel to S(t). Applying the coordinate transformation Eq (38) to an integral over space replaces the infinitesimal volume element dr 0 with where J ¼ det ðdcÞ det @c @u 1 Á Á Á @c @u NÀ1 @c @t ð40Þ is the Jacobian of the transformation Eq (38). The absolute value of this determinant corresponds to the volume of the N-dimensional parallelepiped that has the vectors @c @u 1 ; . . . ; @c @u NÀ1 , and @c @t as adjacent edges. This volume is equal to the volume of the N − 1 dimensional parallelepiped defined by the vectors @c @u 1 ; . . . ; @c @u NÀ1 (base area) multiplied by the projection of @c @t onto the axis perpendicular to this base (i.e., multiplied by the height) [63,65]. Because @c @u 1 ; . . . ; @c @u NÀ1 all belong to the N − 1 dimensional tangent vector space of S(t) at r 0 = ψ(u, t), the unit vector normal to this base is the unit normal vector n, so that the height is j @c @t Á nj. Furthermore, the volume of the N − 1 parallelepiped defined by @c @u 1 ; . . . ; @c @u NÀ1 is equal to is the N × (N − 1) matrix @c @u 1 . . . @c @u NÀ1 [63,65]. With Eq (37), we thus obtain: Since the measure in surface integrals over manifolds is defined as ds ¼ we finally retrieve Eq (35): Note that for the transformation Eq (38) to be injective, it is necessary that its Jacobian is nonzero, and thus that v 6 ¼ 0 in the neighbourhood of r, meaning that no reversal of the direction of propagation of the interface is assumed around r.

Representations of the surface distribution
A similar partition of space Eq (41) can be defined in a neighbourhood of a surface S with N-1 coordinates parallel to S and one coordinate perpendicular to S. Let r 0 ¼ cðu; sÞ 2 R N where ψ(u, 0) is a parameterisation of S with ψ(0, 0) = r 2 S, and with the dependence on s such that @c @s ¼ nðc; sÞ in a small neighbourhood of s = 0. The variable s plays the same role as time t in the developments Eqs (37)- (41), except that it corresponds directly to the arc length along trajectories perpendicular to S, i.e., ds corresponds to |v|dt in Eq (41) and we have dr 0 = dσ ds [18]. This curvilinear partition of space in a small band around S implies in particular that for r 0 in this band: and d S ðr 0 Þ ¼ ðsÞ: Eq (42) represents the factorisation of the Dirac distribution into the coordinates u parrallel to S and the coordinate s perpendicular to S. The denominator accounts for the fact that S is curved. If S is flat and parameterised by orthonormal coordinates, the denominator is one and Eq (42) corresponds (up to a rotation) to the well-known factorisation of the Dirac distribution in cartesian coordinates. It has to be emphasised that for Eqs (42)-(43) to hold, s must be the arc length of a trajectory normal to S. To show Eq (42) we integrate its right hand side over space with a test function φ and use dr 0 = dσ ds: When the surface S is defined implicitly as the zero level of a function φ(r), then 1 O ðrÞ ¼ ðYðrÞÞ, where Θ is the Heaviside step function, and one obtains from Eq (46) the following representation of the surface distribution: This representation of the surface distribution is taken as definition of δ S in [18,20]. It appears in some developments of the level set method [32]. Note that Eq (47) with Eq (43) corresponds to Eq (34) in Section 8.4 of [13].
In [66] the surface distribution appears as the kernel operator where ψ(r) is a parameterisation of S. Indeed, integrating the right hand side of Eq (48) over space with a test function φ gives Z dr Z S dsðuÞ r À cðuÞ ð ÞφðrÞ ¼ One may use also use Eq (42) in Eq (48) to show that it reduces to the representation Eq (43). Eqs (43), (46), (47) and (48) are all different representations of the surface distribution defined by Eq (2). These representations have been used previously in the literature, but were not necessarily identified with a distribution δ S defined by Eq (2) with the properties summarised here, with the notable exception of the early works [18][19][20]. Among the references cited here, Jones [13] probably provides the most rigorous accounts on these representations based on distribution theory, however, without using the suggestive notation δ S(r) .

Balance equation of the indicator function of an evolving domain
The balance equation of the indicator function of an evolving domain can be seen as a particular case of Eq (4) in which η f = 1 and Zðr; tÞ 1 OðtÞ ðrÞ 2 f0; 1g jumps discontinuously between the values 0 and 1, such that: where the regularisation S − (t) on resorption surfaces is implicit in the last two equalities.

Volumetric density of a surface-bound quantity
The surface distribution enables a simple expression for the volumetric density n(r) of a quantity concentrated on a surface S: where ρ is the quantity's surface density on S, and δ S(r) is the surface distribution defined by Eq (2). Eq (51) can be shown by integrating it over a neighbourhood V & R N of a point r on the surface S. The left hand side gives, by definition of n, the absolute amount of the quantity found in the volume V. With the definition Eq (2), the right hand side gives R V\S dσ(p) ρ(p), which by definition of ρ is also the absolute amount of the quantity found on S in the volume V.
Alternatively, the volumetric density of point particles i in space is Assuming the particles all belong to S and using the partition of space Eq (42) and the representation Eq (43), one has nðrÞ ¼ dðsÞ where the surface density on the curved manifold S parameterised by ψ(u) is represented as Appendix B: Numerical discretisation The governing equations for η(r, t) in the one-dimensional and two-dimensional examples were solved numerically based on a simple explicit scheme, using forward finite difference in time (Euler) and a fixed discretisation grid of the computational domain V. The singular surface terms were implemented by explicitly tracking the position of the interface and enforcing the jump condition Eq (3) at this interface. The following steps were performed for each time increment Δt: 1. Evolve the interface given v. Determine the set of discretisation points V f at which η was formed and the set of discretisation points V r at which η was resorbed; 2. For each point r i 2 V f , increase η by η f (r i , t) − η(r i , t); 3. For each point r i 2 V r , set η to 0; 4. For each point r i 2 V, add to η the quantity DtF ðZðr i ; tÞÞ.
In point 2., η(r i , t) is substracted so that the value η f is generated even if η has a residual value at r i . This can happen at reversal points between resorption and formation due to roundoff errors.
In Fig 8, we compare a direct simulation of Eq (18) with the analytical result Eq (20) in the same one-dimensional situation as Fig 3a. The analytic solution Eq (20) requires the arrival timẽ t Ã ðzÞ (in dimensionless coordinate), i.e. the time at which the interface SðtÞ reaches the pointz. In the situation depicted in Fig 8, the arrival time is found numerically by solvingz ¼ Sðt Ã Þ for t Ã using Newton's method, where SðtÞ ¼t þ a sin 2pnt t f with a ¼ 0:35; n ¼ 4;t f ¼ 10.
The double remodelling events simulated in Figs 6 and 7 assumed given populations of osteoblasts and osteoclasts ρ Ob , ρ Oc and constant secretory and resorption rates k f , k r , such that the normal velocity of the interface v = k f Ob > 0 in formation, v = k r Oc < 0 in resorption, was given by: v 2 sin p x À a 2 b 2 À a 2 cos p t À t beg where v 0 = 3 μm/day, t 0 = 25days; v 1 = 0.62 μm/day, t beg to the first remodelling event, and times t beg 2 < t t end 2 to the second remodelling event.