A SIMPLE MODEL OF COLLAGEN REMODELING

. In the present paper we propose and study a simple model of collagen remodeling occurring in latter stage of tendon healing process. The model is an integro-diﬀerential equation describing the possibility of an alignment of collagen ﬁbers in a ﬁnite time. We show that the solutions may either exist globally in time or blow-up in a ﬁnite time depending on initial data. The latter behavior can be related to the healing of injury without the scar formation in a ﬁnite time: a full alignment of collagen ﬁbers. We believe that the present model is an essential ingredient of the full description of collagen remodeling.


1.
Introduction. In the last few decades we have witnessed a rapid growth of use of mathematical methods in the life and medical sciences. The scientific community recognized that the description of complex biological processes with mathematical language is of tremendous importance for their understanding and control, facilitating proper healing, faster cure, regression or at least stabilization of associated diseases. As cancer is one of the major causes of death, particularly in the developed countries, much attention has been paid to mathematical modeling of processes related to cancer disease. The researchers agree that a new field of science-the mathematical oncology-has emerged. However, much less attention has been paid to modeling of processes that do not represent a direct threat to the patient's life. This is the case for traumatic events such as tendon rupture, which, although very frequent and greatly affecting the quality of life, have not been a subject of mathematical modeling. Interestingly, for several decades there was no breakthrough when it comes to treatment of traumatic tendon injuries, to date consisting in a surgery followed by an arduous rehabilitation. Only with the recent developments in therapeutic medicine new prospects for improvement of the tendon healing process have appeared. Still, the use of opportunities offered by therapeutic medicine requires a better understanding of phenomena taking place in healing tendon.
To our knowledge, there is no literature on mathematical modeling of tendon healing processes. Certain similarities or analogies to this process can be found in other types of wound healing already addressed in existing mathematical models. In this respect, literature is quite rich. Early models of wound healing can be found in papers [4,5,18,20]. In 2006 an interesting research on hybrid modeling of fibroblasts migration and collagen deposition during dermal wound healing was published by McDougall at al. [13]. Some review papers on wound healing are for instance [8,19]. However, tendon has its distinctive features, such as very low cellularity and characteristic collagen composition and structure, which renders development of specific tendon-focused models necessary.
In the present paper we consider a simple model of the collagen remodeling occurring in the latter stage of tendon healing process. The model is phenomenological and as such does not describe the underlying biological processes in full detail. Nevertheless, we believe that this is an essential step in proper mathematical description of the collagen fibers remodeling process. We consider an integro-differential equation for which we prove either the global existence or the blow-ups in finite time depending on initial data. The mentioned blow-up may consist in alignment of all collagen fibers to a single orientation, which can be interpreted as injury healing without scar formation (in a finite time). We believe that the model is suitable to describe also other processes where a kind of self-organization occurs.
The paper is organized as follows. First, in Section 2, we give a brief introduction to acquaint the reader with the construction and characteristics of the tendon tissue as well as the specificity of its healing process. In Section 3 we present the model which we next analyze in Section 4. In Section 5 we present the numerical simulations for the proposed model and in Section 6 we provide final conclusions.
2. Biological background. Tendons attaching muscles to bones consists of tightly arranged in parallel collagen (mainly type I) and elastin fibers with less ground substance and relatively low number of tendon specific cells (tenocytes) [14]. Collagen within tendon accounts for 70% of its dry weight and has hierarchical structure of increasing complexity: fibrils, fibers (primary bundles), fascicles (secondary bundles), tertiary bundles and the tendon itself [17]. Special, strongly hydrophilic ground substance surrounding the collagen and cells enables the rapid diffusion of water soluble molecules and cell migration [17]. Tendon cells are arranged in narrow rows laying between thick rows of collagen [14].
The specific structure of the tendon, mainly related to its poor blood supply and low cellularity, causes the healing process following acute injures to be long and fraught with high risk of complications. Nevertheless, the stages of tendon healing are similar to those of other connective tissues, i.e. there are three consecutive, however overlapping phases: i) inflammatory; ii) proliferative or reparative; and iii) remodeling phase [6,17,21].
The inflammation begins immediately after the accident. At the site of the injury blood clot begins to form and starts to activate the release of pro-inflammatory chemicals that, in turn, cause an influx of fluids and inflammatory cells from the surrounding tissue [12]. Inflammatory cells engage in the phagocytosis of devitalized tissue and debris and break down the hemorrhage and clot [3,12,21]. The recruitment and activation of tenocytes begins at that time, however peaks in the proliferative stage [21].
The subsequent proliferative phase begins roughly two days after the injury [21]. Tendon cells continue to migrate and proliferate at the injury site and start to produce collagen fibers of type III that are more randomly oriented then those of type I mainly present in the intact tendon [3,6,21]. The level of some inflammatory cells start to decline. At the same time, other inflammatory cells such as macrophages continue to release growth factors that activate and direct tenocytes recruitment towards the injury site [21].
Eventually, the remodeling phase begins one or two months after the injury [21]. At that time, the type III collagen fibers are replaced with type I collagen fibers that become reoriented and aligned in the direction of stress [21]. From about six till about ten weeks after the injury the reconstructed tissue changes from cellular to fibrous [17]. The remodeling phase of tendon healing may last as long as eighteen months and even surgically repaired tendons never return to their full strength.
3. Model. As we mentioned in the previous section, an essential feature of tendon structure is the collagen fibers orientation, that are aligned in parallel within the healthy tendon. As a result of the tendon injury, particularly its rupture, this parallel structure is disturbed and it is not possible to reconstruct it fully during surgery. The healing process consists mainly in restoration of parallel structure of collagen fibers, however, it is usually possible only to some extent. Scars that form during the healing process fix the improper alignment of collagen fibers. That causes that even after a year the structure and function of the healed tendon remains inferior to the healthy one [21]. Therefore, one of the most important indicators of the successful treatment of tendon injury is the degree of alignment of collagen fibers. In our model, we study the equation describing the simplified process of collagen remodeling that takes place in the latter phase of tendon healing.
The collagen fibers remodeling is a result of action of activated tendon-specific cells, i.e. tenocytes. For simplicity we assume that tenocytes are uniformly distributed within the tendon domain. Therefore we do not consider the density of tenocytes directly, however we would like to point out that in general the density and placement of tenocytes varies during the healing process.
In our model, the function f (t, x, v) describes collagen distribution, i.e. the probability density f = f (t, x, v) to find a collagen fiber at the instant of time t > 0 at point x ∈ D and with orientation v ∈ V, where D, V are bounded domains in R d . We note that the generalization into unbounded domains is straightforward. The variable v, called orientation, gives some information on the alignment of collagen fibers.
We consider the following equation: where T f (y, v; x, w) is the turning rate of orientation w ∈ V at x ∈ D to orientation v ∈ V at x caused by an adaptation to the orientation at y ∈ D. The modeling process leads to the proper choice of the function T f that in general may depend on both collagen distribution f and tenocytes density c. In the present paper however, as already mentioned, we are going to consider a simplified case of constant tenocytes density (uniform distribution) leaving the general case to forthcoming publication. Similar space-independent general structures were proposed by Othmer at al. in [15] and studied e.g. in [7,9] (see also references therein). A detailed bibliographical review can be found in [11]. Various kinetic equations were discussed in [1]. A very simple choice of the turning rate could be That leads to the simple equation Note the "mixed" role of the variable x and v in the equation. The equation seems to be original and never studied in literature. It is somehow similar to those considered in [10,11,16] but its nature is different.
A little more complex, and more realistic, definition where γ > 0, describes the strength of influence of collagen fibers from neighborhood on collagen fiber in considered point. In reality, this influence is mediated by tendon cells, not considered in our model explicitly. The bigger is the γ the stronger the influence is. That leads to the following general class of equations The function β(y, v; x, w) describes the transition from orientation w to orientation v at a point x, as a result of an interaction with another collagen fiber with orientation v located at a point y. From the biological point of view, it would be reasonable to assume that it depends decreasingly on the distance |y − x| and has a small support concentrated close to y = x.
The simplest choice could be γ = 1, however the case γ > 1 leads to more difficult considerations involving possibilities of blow-ups in a finite time, see Theorems 4.3 and 4.4. As we already pointed out, the blow-ups in our model are essential for the description of the tendon healing process.
is Lipschitz continuous, conservative in L 1 (D × V) and its structure is such that the non-negativity of the solutions to Eq. (3) is preserved under the assumption that For any initial datum f 0 being a probability density, for any t > 0, there exists a unique function This result clearly holds for Eq. (2). Consider the following marginal probability densities: By the form of Eq. (2) it is easy to see that the marginal densities corresponding to its solutions are such that and where 0 (x) and ς 0 (v) are the initial distributions Then, Eq. (2) leads to and consequently we have for every t > 0, x ∈ D and v ∈ V. We hence obtain: For any solution f = f (t) to Eq. (2) that is a probability density for t ≥ 0, the marginal probability densities (t, x) and ς(t, v) are conserved, i.e. Eqs. (5) and (6) are satisfied, and where the constant C depends on f 0 .
The latter relation excludes the possibility of blow-ups that represent the successful healing in our model. Therefore, despite its nice mathematical structure, Eq. (2) is not able to describe properly the process in question. The case γ > 1 is much more complex; c.f. simpler equations in [10,11]. In general, we have only local existence . The solution f = f (t) is a probability density and the marginal density is conserved, i.e. (t, for each t ∈ [0, t 0 ]. Moreover, if both β and f 0 are continuous, Proof. The proof is standard. Equation (9) follows by the direct integration of Eq.
(3) with respect to v ∈ V.
The conservation of the marginal density = (t, x), i.e. Eq. (9), is related to the fact that we consider neither movement nor proliferation nor destruction in the model Eq. (3). This property is a priori satisfied even in the general case of Eq. (1).
Assume that β = 1. Let and, as before, ς = ς [1] . Equation (3) may be rewritten as therefore Eq. (9) holds and the marginal density is conserved. Thus we have Introducing By Eq. (11) we have and It follows that ς is not conserved, in contrary to the case covered by Remark 1. The dynamics of ς is governed by the function ς [γ] , hence we analyze its properties. By Eq. (12) we obtain Integration with respect to x yields It follows that Given We have Then, Integration with respect to t yields Denoting (t, The continuity assumption in Theorem 4.3, that certainly may be relaxed, serves to properly define U * .
We are going to provide a more explicit calculation in the case β = 1 and γ = 2. We have then and assume that i.e. ∆(v) < 0 for v ∈ V. Note that ∆(v) ≤ 0, by virtue of the Cauchy-Schwarz inequality. We obtain the explicit expression for H v (U ): Herewith we have the expression for : By Eq. (13), keeping in mind Eq. (14), we have the following ODE: Recall the initial datum ϑ(0) = 0, which results from Eq. (14). Since the righthand side of Eq. (18) is positive, the function ϑ = ϑ(t) is increasing on some interval [0, ϑ * ), where ϑ * := lim t↑t * ϑ(t) can be derived from the limit condition This condition is due to the denominator in the integrand in Eq. (18). This denominator is a continuous function which varies from 1 to −∞, so a singularity appears if it meets zero. We may note that the blow-ups may occur if the denominator in Eq. (18) has a singularity, i.e. Eq. (19) is satisfied. This means that the initial set of orientations is dominated by a "small" number of orientations, giving the proper meaning of the tendon healing process. The situation here is somehow similar to those considered in Refs. [10,11]. Let the right-hand-side of Eq. (18) be denoted by Φ = Φ(ϑ) and Φ(ϑ * ) = lim ϑ↑ϑ * Φ(ϑ) .
We then have-cf. Theorem 2 in [10]-the following result In Items 1 and 2(a), the blow-up time T < ∞ is given by and in Item 2(b) we have T = ∞.

Remark 2.
As already pointed out, contrary to the case covered by Remark 1, ς(t, v) is not conserved even for factorized initial distributions Note that the latter case corresponds to Ξ = 0, unlike in Eq. (16). Then, the solution f (t, x, v) = 0 (x)ς(t, v) leads to an equation similar to the one considered in [10]. Therefore, we may have blow-ups like in [10].
Theorem 4.4 states that for γ = 2, depending on the initial datum, the solutions may be defined on the infinite interval of time (global solutions) or may blow-up in a finite time interval. The latter behavior, as we already mentioned in Section 1, may be related to successful healing.   Simulation results posterior to the blow-up point are not relevant. To handle this issue and filter out irrelevant results, we used a criterion to detect blow-up, which, if triggered, was a sign to terminate the simulation. This criterion was the detection of a high mass concentration. More precisely, the criterion was triggered at the time step t k if there was a node (x m , v n ) of the mesh such that at least 75% of fibers in x m were concentrated in v n . In such cases the simulation was terminated and the previous time step, t k−1 , was assumed to be the last relevant one 1 .
In the first of the simulations, for each x ∈ D, the initial condition f 0 (x, . ) was a combination of two peaks with compact support, each of the peaks having a unique maximum. The height of each peak varied with x independently. The heights and locations of both peaks were such that the resulting initial condition consisted of two lanes, with mass of the one lane dominating the mass of the other, see Figure 1a. Figure 1b and Figure 1c present the evolution of the collagen density and structure in the first simulation. The criterion of 75% mass concentration was triggered at t = 15.4, hence we qualify t = 15.3 as the last relevant step of the simulation, see Figure 1c. For each x ∈ D, the mass of the solution concentrated at a single orientation, probably tending to a blow-up. Since collagen fibers in a healthy tendon are aligned to a single direction, the considered case can be biologically interpreted as a fully successful healing process.
The initial condition in the second simulation was similar, however it is based on combinations of peaks with "truncated tops", hence possessing a level set for the maximal values of positive measure, for each x ∈ D. Thus, unlike the first initial condition, the second one had a plateau, see Figure 2a. Figure 2b and Figure 2c show the evolution of the collagen density and structure in the second simulation. The simulation was performed until the time t = 30.0 and the mass concentration criterion was not triggered within the time span of the simulation. Moreover, the numerical solution achieved a state which underwent no further visible changes, thus probably approximating a stationary state, see Figure 2c. In this state, the mass is distributed uniformly over a span of orientations rather than concentrate in a single orientation. At interpretation level this means that the orientations of fibers are diverse, i.e. the collagen matrix do not regain its natural single-directional structure at full extent. From medical perspective such scenario can be interpreted as formation of the persistent scar which is a frequent complications after a surgery.
The above observations illustrate that, under the dynamics of the model (3) (with β ≡ 1 and γ = 2), the initial arrangement of collagen fibers in the remodeling phase has an essential impact on the final results of the healing process. This highlights the importance of proper surgical treatment, which determines the initial arrangement of fibers. 6. Conclusions. We believe that the existence of blow-ups proven for some initial data and illustrated by numerical simulations can be attributed to the alignment of collagen fibers in a finite time. This can be interpreted as the healing of injury without scar formation. In our opinion, this property of the model is very interesting and deserves much attention as the scar formation is a frequent complication after the tendons injury, greatly affecting their biomechanical properties.
On the mathematical level, we were able to study the simplified case when β is constant. Although it is an essential simplification, this case reflects important properties of the equation. The more realistic case of interaction rate β that has a small support and defines the interaction in a neighborhood is still an open and challenging problem.
The performed mathematical analysis and simulations indicate crucial role of the initial conditions in the final structure of collagen fibers. From the medical perspective, our results stress the importance of the way of performing the intervention after the trauma. The proposed model, although very simplified, turned out to have very interesting properties from both mathematical and modeling perspectives. In the future, we intend to extend the model to include explicitly tendon cells and chemical signaling to allow cells to move and change the structure of collagen fibers. Such extension will lead to the investigation whether the prolonged activation of the tendon cells, achieved through an administration of external stimuli, reduces the occurrence of postoperative scars, therefore leading to better therapeutic outcome.