Stability analysis for a peri-implant osseointegration model

We investigate stability of the solution of a set of partial differential equations, which is used to model a peri-implant osseointegration process. For certain parameter values, the solution has a ‘wave-like’ profile, which appears in the distribution of osteogenic cells, osteoblasts, growth factor and bone matrix. This ‘wave-like’ profile contradicts experimental observations. In our study we investigate the conditions, under which such profile appears in the solution. Those conditions are determined in terms of model parameters, by means of linear stability analysis, carried out at one of the constant solutions of the simplified system. The stability analysis was carried out for the reduced system of PDE’s, of which we prove, that it is equivalent to the original system of equations, with respect to the stability properties of constant solutions. The conclusions, derived from the linear stability analysis, are extended for the case of large perturbations. If the constant solution is unstable, then the solution of the system never converges to this constant solution. The analytical results are validated with finite element simulations. The simulations show, that stability of the constant solution could determine the behavior of the solution of the whole system, if certain initial conditions are considered.


Introduction
A number of models were proposed so far for the process of bone formation. It is reported by many researchers, that mechanical stimulation is an important factor, which influences bone formation. For example, Vandamme et al. (2007a,b,c), Duyck et al. (2007) investigated peri-implant bone ingrowth under well controlled mechanical loading of the interface tissue, and reported that relative implant-interface tissue micromotions qualitatively and quantitatively altered the osseointegration process. The mechanoregulatory models for bone formation were defined, for instance, in Andreykiv (2006), Carter et al. (1998), Claes and Heigele (1999), Doblaré and García-Aznar (2006), Prendergast et al. (1997).
Another biological model for peri-implant osseointegration was proposed in Moreo et al. (2009). This model can be used to simulate osseointegration under a low-medium loading regime taking into account implant surface microtopography. The author did not introduce explicitly the dependence of cell and tissue processes on the mechanical stimulus, and outlined the incorporation of differentiation laws in terms of mechanical variables as one of the future lines of research. The results presented in Moreo et al. (2009) were in agreement with experiments. They predicted that bone formation can occur through contact osteogenesis and distance osteogenesis.
Though, we found that the system of equations, proposed in Moreo et al. (2009), is characterized by appearance of a 'wave-like' profile in the solution for a certain range of parameters. This profile was initially recognized in the solution of the model equations for the 1D domain of length 2.5 mm (Fig. 1b). This domain was chosen for the simulations of bone formation near the cylindrical implant, located within the bone chamber, used in the experiments by Vandamme et al. (2007a,b,c); Duyck et al. (2007). The authors reported, that in a new bone was formed at all distances from the host bone, and integration of bone and implant was achieved. That 'wave-like' profile has not been noticed by Moreo et al. (2009), since for the geometry used in his simulations, in which the distance from host bone to implant was 0.6 mm, only a part of 'wave' is visible in the solution (Fig. 1a), and a 'wave-like' profile is not distinguishable. For larger domains, more 'waves' appear in the solution. The solution for the domain of length 5 mm is shown in Fig. 1c.
The conditions, under which a 'wave-like' profile appears, are studied in the present work. Such a 'wave-like' profile in the solution for cell densities and growth factor concentrations is not realistic. In some cases it also leads to a 'wave-like' distribution of bone matrix inside the peri-implant region. This distribution is in contradiction with experimental observations, from which it follows, that bone forms by deposition on the preexisting bone matrix, and no isolated bone regions appear. Thus, it is desirable to avoid such a profile in the solution of the original model by Moreo et al. (2009), and to take into account the stability properties of the system of equations when introducing mechanical variables in it.
The proposed approach is to study the linear stability of the constant solutions of the system. As the full system of equations is large and extremely complicated for analytic derivations, an equivalent simplified system with similar properties will be defined. The phenomenon of a 'wave-like' profile in the solution could be related to the appearance of bacterial patterns in liquid medium, described mathematically by similar partial differential equations. Those pattern analyses could be found in Myerscough and Murray (1992), Tyson et al. (1999), Miyata and Sasaki (2006).
In Sect. 2 the system of equations proposed in Moreo et al. (2009) is reviewed. The linear stability analysis of the system is carried out in Sect. 3. In Sect. 4 analysis is validated with a sequence of numerical simulations. Finally, in Sect. 5 some conclusions are drawn.

Biological model
The original model proposed in Moreo et al. (2009) consists of eight equations, defined for eight variables, representing densities of platelets c, osteogenic cells m, osteoblasts b, concentrations of two generic growth factor types s 1 and s 2 , and volume fractions of fibrin network v f n , woven bone v w , and lamellar bone v l . The above notations are introduced for non-dimensional cell densities and growth factor concentrations, i.e. for those, related to some characteristic values. Iff and f c are notations of a dimensional variable and of its characteristic value, then a non-dimensional variable f is defined as f =f / f c , f = c, m, b, s 1 , s 2 . The following characteristic values are proposed: c c = 10 8 platelets ml −1 , m c = 10 6 cells ml −1 , b c = 10 6 cells ml −1 , s 1 c = 100 ng ml −1 , s 2 c = 100 ng ml −1 . The model equations are: where D m and A c are the coefficients of random migration and death of platelets. The term ∇ · [H c c∇ p] represents a "linear taxis". It accounts for the transport of platelets towards the gradient of the adsorbed proteins p, which is a predefined function of the distance from the implant surface d. According to Moreo et al. (2009), it is defined as where the terms in the right-hand side represent random migration, chemo-taxis, cell proliferation, differentiation into osteoblasts, and death respectively.
where A b is the rate of osteoblast death.
where the terms in the right-hand side model random migration, growth factor secretion and decay respectively.
where the first term in the right-hand side determines random migration, the second and the third ones-growth factor secretion, and the last one-decay.
Remark 1 Growth factor 1 s 1 is assumed to stimulate the differentiation of osteogenic cells into osteoblasts. In Moreo et al. (2009) originally, the differentiation term in Eqs.
(2) and (3) was given in the form α mb s 1 β mb +s 1 m. In this paper, we propose a more general representation for differentiation, which is (α p0 + α mb s 1 β mb +s 1 )m. Parameter α p0 implies, that differentiation can take place, if the growth factor 1 concentration s 1 is zero. This assumption is not in contradiction with experimental observations (Linkhart et al. 1996;Dimitriou et al. 2005), and can be useful, in order to get a more realistic simulation results for different problems. The profit of this representation for differentiation will be demonstrated in Remark 3 in Sect. 3.1.
We also propose the following values, assuming, that differentiation also takes place without presence of the growth factor: In our present study, we will consider both sets of parameter values in Eqs. (12) and (13).

The simplified biological model
Our present aim is to study the conditions, that give the appearance of a wave-like profile. Simulations, performed for the full system, show that the wave-like profile can appear in the solution for densities of osteogenic cells m and osteoblasts b, for growth factor 2 concentration s 2 , and for volume fractions of fibrin network v f n , woven bone v w and lamellar bone v l , if the computational domain is sufficiently large. The equations for variables m, b and s 2 (2), (5), (3) are coupled and can be solved, after the solution for c and s 1 is obtained from the Eqs. (1) and (4). The equations for variables v f n , v w and v l (6), (7), (8) contain only reaction terms in their right-hand side. The wave-like profile in the solution for these variables appears due to the wave-like profile in the solution for osteoblasts and growth factor 2. Therefore we will study the phenomenon of the wave-like profile in the solution for variables m, b and s 2 . The solution for m, b and s 2 is determined by the system of Eqs. (1)-(5).
We assume, that the profile appearance could be related to the stability of the constant solutions of the system. Zero solutions c = 0, s 1 = 0 are the only constant solutions of system (1)-(5) for variables c and s 1 .
The equations for platelets c and growth factor 1 s 1 (1) and (4), can be solved separately from the other equations. That means, that the evolution of the platelet density c(x, t) and growth factor 1 concentration s 1 (x, t) does not depend on the evolution of other biological and chemical species involved in the model. Equation (1) contains a term, corresponding to the death of platelets, but it does not contain a term, corresponding to the production of platelets. Therefore, the total amount of platelets decays to zero with time. The production of growth factor 1 s 1 is proportional to platelets concentration, and thus the production of s 1 also decays with time, while death rate A s1 is constant in time. It can be proved, that the integrals of platelet density and growth factor 1 concentration over the problem domain tend to zero with time, if zero flux on the boundaries is considered. If negative values in the solution for c(x, t) and s 1 (x, t) are avoided (otherwise the solution becomes biologically irrelevant), then it follows, that these functions tend to zero almost everywhere in the problem domain. Numerical simulations confirm (Fig. 2), that for a large time t the solution s 1 (x, t) is very close to zero. The stability analysis deals with the asymptotic behavior of the system, that is with the behavior of the solution for long time periods. Therefore, we derive the simplified system from Eqs. (2), (5) and (3), assuming s 1 (x, t) ≡ 0, which gives Remark 2 In deriving (15) we assumed, that α b2 = α m2 and β b2 = β m2 . These simplifying assumptions are in line with the values for α b2 , α m2 , β b2 and β m2 , proposed by Moreo et al. (2009), which were introduced in (11).
Remark 3 As it was mentioned, the concentration of s 1 becomes close to zero after a certain period of time. Then, differentiation of osteogenic cells into osteoblasts is roughly described by the term α p0 m, as this is done in Eqs. (14), (16). This term turns to zero, if α p0 = 0, as was proposed by Moreo et al. (2009). Solution of (16), defined as b(x, t) = b 0 (x)e −A b t , converges to zero with time. From biological point of view, this means, that osteogenic cells stop to differentiate after a certain time period. There is no source of newly formed osteoblasts, and their amount decrease to zero, due to cell death. If α p0 = 0, then differentiation takes place also when s 1 is zero. This allows to obtain the solution for osteoblasts, which does not converge to zero, and hence, is more realistic from biological point of view. For this reason, we consider both the parameter values in Eq. (12), as proposed by Moreo et al. (2009), and the alternative values in Eq. (13). Moreo et al. (2009) investigated the linear stability of the constant solutions of the system, which is similar to system (14)-(16), against purely temporal perturbations. In this paper we will study the system stability against arbitrary perturbations (also non-homogeneous perturbations).
We mention here, that for the existence of real s 2± the necessary condition is: This necessary condition is written in terms of the model parameters as: From (24) it is derived, that (23) is equivalent to: The sign of s 2± depends on the sign of coefficients a 1 and a 0 (coefficient a 2 is larger than zero, which follows from its definition). Both roots will be positive if a 1 < 0 and a 0 > 0 and if inequality (23) holds.
For the parameter values in Eqs. (11), (12) the constant solutions have values: Remark 5 For the chosen parameter sets (11), (12) and (11), (13), growth factor 2 concentration s 2− is negative, which is unphysical. It is desirable to avoid such a negative concentration of growth factor 2 in the solution of the problem (14)-(16). Calculations show, that for the chosen parameter values there are two positive eigenvalues of the Jacobean of the equation system, linearized for the case of small purely temporal perturbations near the constant solution z − . Hence, constant solution z − is unstable against temporal perturbations. In simulations we were able to avoid negative values in the solution for s 2 , by choosing sufficiently small time step and mesh size and starting with positive initial values for concentrations of cells and growth factor.

Non-homogeneous perturbations
Next we propose an approach, to study the stability of constant solutions of the system (14)-(16). This approach is valid for a domain in any coordinate system, for which eigenfunctions of Laplace operator can be found. In this paper, the examples of the eigenfunctions are given for the domains in 1D Cartesian coordinates and in axisymmetric coordinates, which have one independent coordinate. The independent space coordinate is denoted by x for both coordinate systems. Suppose that non-homo- Then the solution is given in the form: where |ε| 1. Then we substitute (27) into (14)- (16), and linearize with respect to small ε: Let us denote the problem domain as [x 0 , x 0 + L]. Assume, that on the boundaries the flux of cells and growth factor is zero. Then we consider perturbations of the form: and considered boundary conditions, i.e. zero flux on the boundaries: ∇φ n (x 0 ) = ∇φ n (x 0 + L) = 0.
If Cartesian coordinates are considered, then the function φ n (x) is given as φ C n (x) = cos(k n (x − x 0 )), where k n = π n L , n = 1, 2, . . .. In this case k n is a wavenumber.
In the case of axisymmetric coordinates functions φ n (x) have the form φ a Functions φ a n (x), n = 1, 2, . . . are not periodic. They could be roughly described as "waves" with variable in space wavelength and magnitude. For simplicity, k n will be referred to as 'wavenumber', also if it is introduced in functions φ a n (x).
Substituting (29) into (28), we get: where where Then from (30): where C 0 n define the perturbations imposed on the constant solution of the system initially at time t = 0: Thus the solution of (28) is written as: The magnitude of perturbations C n (t) = e A kn t C 0 n of mode n, will grow in time, if at least one of the eigenvalues of matrix A k n is a positive real number or a complex number with a positive real part. And C n (t) will converge to zero, if all the eigenvalues of A k n are real negative, or complex numbers with the real part less than zero. If the matrix A k n has precisely one zero eigenvalue, and other eigenvalues are real negative of complex with negative real part, then small perturbations remain small for infinite time period.
It is not complicated to find expressions for the eigenvalues of A k n , evaluated at the 'chronic non healing state' z t = (0, 0, 0) and 'low density state' z 0 = (m 0 , 0, b 0 ). For the constant solution z t eigenvalues of A k n are: Therefore, if m 0 is positive, constant solution z t is unstable against purely temporal perturbations and perturbations with small wavenumber 0 < k n < α m0 m 0 D m . The first eigenvalue λ 1t (k 2 n ) takes the largest positive value for wavenumber k 0 , i.e. for the purely temporal perturbation mode.
Remark 7 If we consider negative m 0 , then 'chronic non-healing state' z t will become stable against perturbations with any wavenumber. Further the constant solution z 0 will contain an unphysical negative concentration for osteogenic cells.
< 0 implies, that differentiation and death of osteogenic cell dominate over their production. Therefore, this situation is not relevant for the considered model of bone formation, and further m 0 > 0 is assumed a priori.
For the constant solution z 0 = (m 0 , 0, b 0 ) matrix A k n eigenvalues are: takes a positive value, which is true for the current parameter values in Eqs. (11), (12) and (13), then the constant solution z 0 is unstable against perturbations with wavenumbers k 2 The largest eigenvalue λ 20 corresponds to zero wavenumber k 0 , i.e. to the purely temporal mode of perturbation.
The eigenvalues of matrix A k n defined at points z − and z + could not be found in such a trivial manner, as for constant solutions z t and z 0 . They are obtained from the characteristic equation, which is a non-trivial cubic algebraic equation. Therefore, instead of analyzing the expressions for the eigenvalues, which are extremely complicated in this case, we propose a different approach, based on a reduction to two equations with similar stability properties, to study the stability of the considered system of equations.
Remark 8 For the chosen parameter values, see expressions (11), (12) and (11), (13), s 2− is negative, hence constant solution z − is biologically irrelevant in this case. Therefore, we will only analyze the stability of constant solution z + and not of z − . The stability analysis, being introduced for z + , is not valid for the constant solution z − , if it contains the negative value of growth factor concentration. Calculations also show, that for parameter values (11), (12) and (13), constant solution z − is unstable against at least purely temporal perturbations.

Stability of the system of two equations
To simplify the stability analysis, system (14)-(16) is reduced to a system of two equations. For this reduced system the assumption is made, instead of Eq. (16). The system is defined as: Substitution of (37) into Eq. (16), yields the condition ∂b ∂t = 0, which is not true in general case. Therefore, system (14)-(16) and system (38) are not equivalent, and their stability properties are different in general. However, it will be shown in Sect. 3.4, that there is a certain similarity (or correspondence) between the stability properties of the two systems. This similarity is sufficient, to transfer important results, obtained from the stability analysis for the system of two Eqs. (38), onto the system of three Eqs. (14)-(16).
System (38) has constant solutions, that are analogous to those of the system (14)-(16). They are: Considering solutions of the form and substituting them into system (39), for each n = 0, 1, . . . we arrive at: First we investigate the stability properties of system (39) and then determine, how they are related to the stability properties of the system of three Eq. (28). Since s 2+ = −β m2 , then from (20) m + = 0. Therefore, matrix A k n , evaluated at point (m + , s 2+ ), can be simplified. From the first equation of the system (17) we get: Then: where χ is defined in (22). Considering (20), we arrive at Everywhere in the calculations, presented in Moreo et al. (2009) and in this paper, the same values are used for parameters β m and β m2 . So both notations β m and β m2 is used, though β m2 = β m is supposed below. Then Therefore, we end up with Then the characteristic equation for matrix A k n , evaluated at point (m + , s 2+ ), is given by: From Eq. (41) the eigenvalues of A k n (m + , s 2+ ) are determined as: We mention that, if Thus, we can formulate the lemma. The wavenumbers which lead to growing perturbations are determined by inequality c(k 2 n ) < 0. We can write c(k 2 n ) is the form: where  (19), then γ 0 defined in (47) is non-negative.
We mention here, that under assumptions of Lemma 2, c(0) = γ 0 ≥ 0. Then from Lemma 1 we arrive at Lemma 3. (18) is positive, β m2 = β m and there exists a real positive s 2+ defined in (19), then for zero wavenumber k 0 , matrix A k n (m + , s 2+ ) has either one zero eigenvalue and one negative, or two negative eigenvalues, or two complex eigenvalues with negative real part; and the constant solution (m + , s 2+ ) of the system (38) is stable against the purely temporal perturbations.

Lemma 3 If, for the chosen parameter values, m 0 defined in
Since k n ∈ [0, ∞), then c(k 2 n ), given in (44) could be considered as a real function of a real non-negative argument. It is a quadratic polynomial. The interval, where c(k 2 n ) < 0, is defined by the roots of the polynomial. If this polynomial has no roots among non-negative real numbers, then for ∀k n ∈ [0, ∞), c(k 2 n ) > 0, since γ 2 defined (45) is positive. Thus, it is necessary to find the conditions, if polynomial defined in (44) has at least one non-negative real root. The general formula for the roots of the polynomial is: The discriminant of the polynomial is: Since γ 2 > 0 and γ 0 ≥ 0 under the conditions of Lemma 2, the polynomial c(k 2 n ) has either two real roots of the same sign as −γ 1 , which are different if D γ > 0, and coincident if D γ = 0; or two complex roots with real part − γ 1 2γ 2 , if D γ < 0. (18) is positive, β m = β m2 and there exists a real positive s 2+ defined in (19). Then if D γ defined in (55) is positive, and γ 1 defined in (54) is negative, then ∃κ 1 , κ 2 ∈ R defined by expression (52), such that 0 ≤ κ 1 < κ 2 , and the constant solutionz + = (m + , s 2+ ) of the system (38) is unstable with respect to the perturbations with wavenumbers k n ∈ (κ 1 , κ 2 ). Otherwise, constant solutionz + is stable.
Otherwise, ∀k n ∈ [0, ∞) the eigenvalues of matrix A k n are either real non-positive numbers (the matrix A k n can not have more than one zero eigenvalue) or complex numbers with negative real part. Hence, initially small perturbations remain small during any period of time, or even disappear when t → ∞, and constant solutionz + is stable in this case.
The parameters γ 1 and D γ , can be written in terms of model parameters as

Correspondence between the systems of two and three equations
Next we determine the relation between the eigenvalues of matrices A k n (m + , s 2+ ) and A k n (m + , s 2+ , b + ), in order to demonstrate the similarity between the stability of systems (14)- (16) and (38), with respect to perturbations about equilibria (m + , s 2+ , b + ) and (m + , s 2+ ) respectively. Let us define matrix M k n : From the definition of A k n we have: A k n (2,3) = A k n (2,1) . Then The determinant of this matrix is the characteristic polynomial of A k n : From the definition of matrices A k n and A k n , it follows that A k n (1,1) = A k n (1,1) , (2,2) . Therefore, the determinant of matrix A k n − λI 2 and characteristic polynomial of matrix A k n is From (56) and (57) we derive: 2,1) . (58) Then we denote the characteristic polynomials of matrices A k n and A k n , which are evaluated at the constant solutions (m + , s 2+ , b + ) and (m + , s 2+ ) respectively, as cubic polynomial P 3 (λ) and quadratic polynomial P 2 (λ) with respect to λ: Equation (58) could be written as: where If s 2+ > 0, it follows from (20) that m + > 0, and from (40) we obtain: Lemma 4 Suppose, that for the chosen parameter values m 0 defined in (18) is positive, and that there exists a real positive s 2+ defined in (19). If the matrix A k n (m + , s 2+ ) has one real negative eigenvalueλ 1 < 0 and one real positive eigenvalueλ 2 > 0, then A k n (m + , s 2+ , b + ) has one real positive eigenvalue and either two real negative eigenvalues, or two complex conjugated eigenvalues with negative real part.
Proof From the assumption of the lemma and from (61) it follows, that C(k 2 n ) > 0. Let A k n (m + , s 2+ ) have one real negative eigenvalueλ 1 < 0 and one real positive eigenvalueλ 2 > 0. The characteristic polynomial can be written as P 2 (λ) = (λ −λ 1 )(λ −λ 2 ). Then from Eq. (59) From (62) we get: Since P 3 (λ) is continuous, it follows from (63), that polynomial P 3 (λ) has at least one real positive root λ 1 on the interval (0,λ 2 ). The other two eigenvalues λ 2 and λ 3 of A k n (m + , s 2+ , b + ) could be real (negative or positive) or complex conjugated numbers (as the coefficients of the polynomial are real). We can write: since this polynomial has λ 1 , λ 2 , λ 3 as its roots. As the coefficients at the second degree of λ in the two expressions for P 3 (λ) from (62) and (64) should be equal, we have λ 2 + λ 3 =λ 1 +λ 2 − A b − λ 1 . From (42) it is derived: The above inequality holds, since it was mentioned in (43), that b(k 2 n ) > 0, if m 0 > 0 and s + > 0. Thus, if two other eigenvalues are real, then from (65) it follows, that at least one of them is negative. Let us suppose λ 2 < 0. Then lim λ→−∞ P 3 (λ) = ∞, and P 3 (0) = −λ 1λ2 A b > 0. That means that on the interval (−∞, 0) polynomial P 3 (λ) does not change its sign, or changes it twice. Since P 3 (λ) is continuous, it follows from λ 2 < 0 that λ 3 also is negative. In the case, when λ 2 and λ 3 are complex conjugated, their real part is λ re = (λ 2 + λ 3 )/2 < 0. (19). If A k n (m + , s 2+ ) has one zero eigenvalue and one real negative eigenvalue, then A k n (m + , s 2+ , b + ) has one zero eigenvalue and either two real negative eigenvalues, or two complex conjugated eigenvalues with negative real part.

Lemma 5 Suppose, that for the chosen parameter values there exists a real positive s 2+ defined in
Proof From the assumption of the lemma and from (61) it follows, that C(k 2 n ) > 0. Let A k n (m + , s 2+ ) have one zero eigenvalue and one real negative eigenvalue, λ 1 <λ 2 = 0. Then the characteristic polynomial P 2 (λ) has the form P 2 (λ) = λ(λ − λ 1 ). Then Eq. (59) implies And eigenvalues of A k n (m + , s 2+ , b + ) are following: Since C(k 2 n ) −λ 1 A b > 0 and A b −λ 1 > 0, then from (67) it follows, that eigenvalues λ 2,3 are either real and negative (possible coincident), or complex with negative real part. (19). If A k n (m + , s 2+ ) has two real negative eigenvalues, then A k n (m + , s 2+ , b + ) has either three real negative eigenvalues, or one real negative eigenvalue, and two complex eigenvalues with negative real part.

Lemma 6 Suppose, that for the chosen parameter values there exists a real positive s 2+ defined in
Proof From the assumption of the lemma and from (61) it follows, that C(k 2 n ) > 0. Let A k n (m + , s 2+ ) have two real negative eigenvaluesλ 1 ≤λ 2 < 0. Then the characteristic polynomial P 2 (λ) has the form P 2 (λ) = (λ −λ 1 )(λ −λ 2 ). Then from Eq. (59) From (68) we get: Since P 3 (λ) is continuous, it follows from (69), that polynomial P 3 (λ) has at least one root on the interval (−A b , 0). Thus we can suppose, that −A b < λ 1 < 0. From (68) it follows, that for λ ≥ 0 polynomial P 3 (λ) only takes values less than zero. That means, that P 3 (λ) has no non-negative real roots P 3 (λ). Thus, if two other eigenvalues of A k n (m + , s 2+ , b + ) are real, they are also negative. Though it is possible, that polynomial P 3 (λ) has two complex conjugated roots. Let us denote them as λ 2,3 = λ re ±iλ im . Then: As the coefficients at the second degree of λ in two expressions for P 3 (λ) (68) and (70) should be equal, we derive: 2λ re =λ 1 +λ 2 − A b − λ 1 . Asλ 1 ≤λ 2 < 0 and −A b − λ 1 < 0, we get that λ re < 0. That is, if two eigenvalues of A k n (m + , s 2+ , b + ) are complex, then their real part is less than zero. (19). If A k n (m + , s 2+ ) has two complex conjugated eigenvalues with negative real part, then A k n (m + , s 2+ , b + ) has either three real negative eigenvalues, or one real negative eigenvalue, and two complex eigenvalues with negative real part.

Lemma 7 Suppose, that for the chosen parameter values there exists a real positive s 2+ defined in
Proof From the assumption of the lemma and from (61) it follows, that C(k 2 n ) > 0. Let A k n (m + , s 2+ ) have the complex conjugated eigenvalues with negative real part λ 1,2 =λ re ± iλ im ,λ re < 0. Then the characteristic polynomial P 2 (λ) takes positive values for ∀λ ∈ R and has the form P 2 (λ) = (λ 2 − 2λ re λ +λ 2 re +λ 2 im ). Then from Eq. (59), we obtain From (71) we get: Since P 3 (λ) is continuous, it follows from (72), that polynomial P 3 (λ) has at least one root on the interval (−A b , 0). Thus we can suppose −A b < λ 1 < 0. From (71) it follows, that for λ ≥ 0 polynomial P 3 (λ) takes values less than zero. That means, that P 3 (λ) has no non-negative real roots P 3 (λ). Therefore, if the two other roots of P 3 (λ) are real, they are also negative.
Next we investigate the possibility, that polynomial P 3 (λ) has two complex conjugated roots. We denote them as λ 2,3 = λ re ± iλ im . Then: As the coefficients of λ 2 in two expressions for P 3 (λ) (71) and (73) should be equal, we derive: 2λ re = 2λ re − A b − λ 1 . Asλ re < 0 and −A b − λ 1 < 0, we get that λ re < 0. That is, if two eigenvalues of A k n (m + , s 2+ , b + ) are complex, then their real part is less than zero.
3.5 Stability of the system of three equations Proof From Lemma 3, 5, 6 and 7, we obtain, that for zero wavenumber k 0 , matrix A k n , evaluated at the constant solution z + , has either two negative eigenvalues and one zero eigenvalue, or three real negative eigenvalues, or one real non-positive eigenvalue, and two complex eigenvalues with negative real part. (18) is positive, β m = β m2 and there exists a real positive s 2+ defined in (19). Then if D γ defined in (55) is positive, and γ 1 defined in (54) is negative, then ∃κ 1 , κ 2 ∈ R defined by expression (52), such that 0 ≤ κ 1 < κ 2 , and the constant solution z + of the system (14)-(16) is unstable with respect to the perturbations with wavenumbers k n ∈ (κ 1 , κ 2 ). Otherwise, constant solution z + is stable.

Theorem 2 Suppose, that for the chosen parameter values m 0 defined in
Proof The theorem is proved, by analogy with the proof of Theorem 1. In that proof all possible cases for the signs of the parameters D γ and γ 1 are considered, and the relations between the eigenvalues of matrix A k n and wavenumber k n are determined for each case. From those relations, and from relations between eigenvalues of matrices A k n and A k n , stated in Lemma 4-7, it is possible to determine the correspondence between the eigenvalues of matrix A k n and wavenumber k n for the sets of signs of the parameters D γ and γ 1 . Therefore, it is obtained, that if D γ > 0, and γ 1 < 0, then ∃κ 1 , κ 2 ∈ R defined by expression (52), such that 0 ≤ κ 1 < κ 2 , and the magnitude of the perturbation modes with wavenumbers k n ∈ (κ 1 , κ 2 ) grow monotonically after a certain period of time, since one of the eigenvalues of A k n is positive. Hence, constant solution z + = (m + , s 2+ , b + ) is unstable with respect to these perturbation modes.
Otherwise, ∀k n ∈ [0, ∞) the eigenvalues of matrix A k n are either real non-positive numbers (the matrix A k n can not have more than one zero eigenvalue) or complex numbers with negative real part. Hence, initially small perturbations remain small during any period of time, or even disappear when t → ∞, and constant solution z + is stable in this case.
The conditions on parameters, stated in Theorem 2, can be formulated in compact form: From the proof of the theorem, it follows, that condition (74) is a necessary condition for the instability of the solution z + , since it is equivalent to the existence of real positive numbers κ 1 and κ 2 . The necessary and sufficient condition holds, if there exist wavenumbers k n ∈ (κ 1 , κ 2 ). From (52) it follows, that the length of interval (κ 1 , κ 2 ) is equal to D γ D m D s2 . If D γ is small enough, then it is possible, that no wavenumber k n will lie inside interval (κ 1 , κ 2 ), and perturbations will not grow. In this case, the necessary condition for instability holds, but the solution is stable.
The necessary instability condition (74), can be transformed into the sufficient stability condition by the substitution of the sign in inequality (74) by the opposite one: This condition is formulated in terms of model parameters and does not depend on the problem statement. This means, that general instruction on the choice of the parameter values, which guaranty the stability of the constant solution z + , can be formulated.
The necessary and sufficient stability condition is opposite to the necessary and sufficient instability condition, which depend on the wavenumbers k n . Set of wavenumbers k n contains infinite number of elements, and is determined by the domain size, by the coordinate system and by the boundary conditions. Therefore, it is not possible to state the necessary and sufficient condition in terms of model parameters for the general case. For different boundary conditions, coordinate systems or domain sizes, these conditions have to be reformulated.

Parameter choice and stability
In this subsection the choice of parameter values, providing stability of the constant solution z + = (m + , s 2+ , b + ) of the system of three equations, is discussed.
For the parameter values in Eq. (11), (12) and (13), the constant solutions z t = (0, 0, 0), z 0 = (m 0 , 0, b 0 ) and z − = (m − , s 2− , b − ) are unstable, and the solution will not converge to these constant solutions. From a biological point of view, this is a favorable situation. Since, the 'non-healing state' z t contains zero concentrations of osteogenic cells and osteoblasts, the 'low density state' z 0 corresponds to much lower concentrations of osteogenic cells and osteoblasts, compared to those for z + , and the constant solution z − contains unphysical negative concentrations of cells.
For the chosen parameter value sets in Eqs. (11), (12) and (11) (12) and (13), and for any B m2 and D m , z − is unstable against purely temporal perturbations. Therefore, varying B m2 and D m , we can achieve the stability of constant solution z + , while constant solutions z t , z 0 and z − remains unstable; -calculations showed, that stability condition (75) is most sensitive with respect to the parameters B m2 and D m . That is, the ratio of the initial parameter value and the ultimate value of the parameter, which satisfies condition (75), is much smaller for B m2 and D m , compared to the rest of the model parameters.
The first quadrant of plane (D m , B m2 ), which contains all possible non-negative values D m and B m2 , can be divided into three regions, with regard to the stability of the solution z + : region R 1 : sufficient stability condition (75) holds, solution z + is stable; region R 2 : condition (75) does not hold, no wavenumbers k n lie in the interval (κ 1 , κ 2 ), solution z + is stable;  (D m , B m2 ), where constant solution z + is stable (R 1 , R 2 ) and unstable (R 3 ), for the case of zero flux of m and s 2 on the boundaries, 1D Cartesian coordinates, and domain of length 0.6 mm. The rest of model parameters are initialized: a as in (11), (12), and b as in (11), (13) region R 3 : condition (75) does not hold, some of wavenumbers k n lie in the interval (κ 1 , κ 2 ), solution z + is unstable.
Configuration of regions R 2 and R 3 depend on the specified boundary conditions, coordinate system and domain length. In Fig. 3 these regions were plotted for the case of zero flux of m and s 2 on the boundaries, 1D Cartesian coordinates, and domain of length 0.6 mm. This length is equal to the width of the domain, used in the numerical simulations by Moreo et al. (2009). Values of model parameters, given in (11), (12) (Fig. 3a), and in (11), (13) (Fig. 3b), were chosen.
With use of (45), (46) and (47), sufficient stability condition (75) can be rewritten as follows: or in the form , and γ 0 is defined in (47). Inequalities (76) and (77) determine the values of B m2 and D m , which ensure the stability of the solution z + .
The following remark can be helpful for the solution of practical problems. Suppose, that initial values of model parameters do not satisfy sufficient condition (75) for the stability of the solution z + . Then, it is possible to guaranty the stability of z + in general case (i.e. for any set of wavenumbers, which are determined by problem statement), by decreasing the value of B m2 , or increasing D m , until condition (76) or condition (77) is satisfied, respectively.

Numerical results
The predictions from the linear stability analysis are validated against a sequence of numerical simulations. The sufficient stability condition is considered in the form (76) and parameter B m2 is varied.
If we fix the values of all parameters, except B m2 , then the right part of inequality (76) can be denoted as the ultimate value B lim m2 , such that for B m2 ≤ B lim m2 small perturbations near (m + , s 2+ , b + ) are predicted not to grow with time. For B m2 > B lim m2 small perturbations of mode φ n (x) will grow, if κ 1 < k n < κ 2 . If B m2 is close to ultimate value B lim m2 , then the interval (κ 1 , κ 2 ) is small, and it could happen, that no wavenumber k n lies inside this interval. In this case perturbations near the constant solution will not grow, in spite of the fact, that sufficient stability condition (76) does not hold.
In Fig. 4 the results of numerical simulations are shown. The solutions were obtained with use of the finite element method. Linear 1D elements of size 0.02 mm were used for the discretization in space. The implicit backward Euler method, to prevent instabilities due to numerical time integration, and adaptive time stepping were used for time integration. Zero flux of m, s2 on the boundaries was specified as the boundary conditions. To introduce the perturbations in the initial solution during simulations, the corresponding constant solution value plus a small random number were assigned to every degree of freedom at time t = 0. From Fig. 4 it follows, that for values B m2 less than the ultimate value, the numerical solution tends to the constant solution (m + , s 2+ , b + ) with time (Fig. 4a). And if parameter B m2 is larger than B lim m2 and such, that ∃k n ∈ (κ 1 , κ 2 ), then there is no convergence to the constant solution, and a wavelike profile occurs in the solution (Fig. 4c, d). However, if B m2 is larger than B lim m2 , but such that no wave number lies inside (κ 1 , κ 2 ) yet, then the numerical solutions again converge to the constant solution (m + , s 2+ , b + ) (Fig. 4b). Thus, the predictions of the linear stability analysis are fully confirmed by the numerical simulations.
The introduced linear stability analysis allows to assess the stability of the considered constant solution. From its stability it can be concluded, whether or not small perturbations grow with time. The important conclusion can be made, for cases in which perturbations are large: if the constant solution is not stable, then the solution of the problem will never converge to that constant solution. Hence, the introduced   (11), (12) linear analysis provides important results also for the case of large perturbations, since it allows to avoid the situation, in which the solution, which is constant in time and in space, can never be reached. However, if the constant solution is stable, it is still unknown, how large initial perturbations behave, whether they disappear or prevail, or even grow.
The simplified system (14)-(16), and the full system (1)-(8) were solved numerically for initial and boundary conditions (9), (10) and (78), (79) respectively, and for a number of parameter value sets. Some of the solutions for the full system (1) (11), (13) after a certain period of time (Fig. 5a). Though, if the constant solution (m + , s 2+ , b + ) is not stable, then a wave-like profile develops in the solution for osteogenic cells and growth factor 2 and for parameter values (11), (13) also in the solution for osteoblasts.
For some values of parameter B m2 that 'wave-like' profile is steady (Fig. 5b). Though, if B m2 is much larger than the ultimate value, then the waves in the numerical solution are not steady, but moving (Fig. 5c, d). This is in agreement with the stability analysis presented in Sect. 4.

Conclusions
We have defined a simplified system of three equations, characterized by the appearance of a wave-like profile in the solution under the same conditions, as for the solution of the full system of eight equations. For the considered parameter values the simplified system has four constant solutions. The sufficient stability condition for one of the constant solutions, denoted as z + = (m + , s 2+ , b + ), is derived in terms of model parameters, by means of the linear stability analysis. If all constant solutions are unstable, then by changing the values of the model parameters B m2 and D m , it is possible to make the solution z + stable, while three other constant solutions z t , z 0 and In real simulations for the peri-implant osseointegration, initial conditions correspond to the large deviations from the constant solution. However, linear stability analysis provides important results also in this case. It allows to avoid such values of model parameters, for which all constant solutions are unstable, and consequently, can not be reached. Linear stability analysis makes it possible to determine parameter values, for which the solution of the problem will never converge to the solution, which is constant in time and in space. This conclusion is confirmed by the numerical simulations, which evidence, that a wave-like profile appears in the solution, if all the constant solutions are unstable. The numerical simulations also show, that if the solution z + is stable and z t , z 0 , z − are unstable, then numerical solutions for unknowns m(x, t), s 2 (x, t), b(x, t) of the full and the simplified system converge to the constant solution (m + , s 2+ , b + ) after a certain period of time, when starting with initial conditions proposed in Moreo et al. (2009). Therefore, the numerical simulations demonstrate, that if constant solutions z t , z 0 , z − are unstable, then stability of the constant solution z + could determine the behavior of the solution of the whole system, if specific initial conditions are considered. That makes it possible to determine the values of model parameters, for which biologically irrelevant solutions with a 'wave-like' profile can be obtained.