Bone remodeling: Komarova model

Dear Community,

I make a code public for bone remodeling, using the Komarova model. It was developed by me alone. I provide the code as is, without warranty.

The agreement between my data and the ones from a reference are very good for fixed loading, but deviate for a more intricate loading. Lacking access to Abaqus with a Fortran compiler, I am at present unable to pinpoint the cause for the discrepancy between my data and data from the reference. The code can serve as a starting point and I give directions how to identify the cause for that discrepancy.

The self-coded part of the subroutine ends at the subroutine ODEX. This ODEX-code provides a numerical solution of a system of first order ordinary differential equations and was taken off the shelf. That part of the program was coded by mathematicians. Don't let it bother you, it acts as black box. Mathematical advisory needed !

I recommend to first study the reference [1] as it presents reference solutions for fixed stimulus values.

List of State Variables:
01 number of osteoclasts
02 number of osteoblasts
03 number of inactive osteoclasts
04 number of inactive osteoblasts
05 number of active osteoclasts
06 number of active osteoblasts
07 density
08 temporal variation of the density
09 variation of the density
10 Young's modulus
11 Poisson's ratio
12 density of the elastic energy
13 stimulus DELTA_E
14 g12
15 g21
16 speed of the variation of the number of osteoclasts = d(n1)/dt
17 speed of the variation of the number of osteoblasts = d(n2)/dt
18 omega (defined in [2] after Eq. (11))
19 zeta (defined in [2], Eq. (13))

 
 
The following simulations are included in the attachment:

1) the simulations V01, V02, V03. These are the three simulations to [1] with fixed values for the stimulus, defined in the subroutine. Thus, loading or deformation of the FEM model was not required and the input-files to V01, V02, V03 are identical.
The Fortran-files differ only in their definition of DELTA_E.

2) V07 defines DELTA_E as given in [2], equation 18. The simulation was not analyzed as V08, with results easier to verify, generated a discrepancy to the expected values, see below. Also, the loading of the model is not as indicated in [2].

3) V08 imposes the stimulus DELTA_E in Fig. 7 from [1]. This stimulus, of course, should follow from external loading of the model. As a preliminary version, the stimulus is defined here as a function of time within the subroutine. The figure referred to was digitized, and a time function allows to interpolate linearly between the time points. Thus, once again, no external loading or deformation of the model was required.

The initial Young's modulus in this simulation follows from the fit formula for the modulus vs. density with the indicated value for the density. There is obviously a misprint in my input-file, the value is presumably a minor error.

The density of an element in my simulation V08 deviates from the one depicted in Fig. 7. This deviation might originate from an inaccuracy in my code. Another possibility is that the authors inadvertently inserted a graph for density vs. time from another simulation. Mind that the strain for that simulation is given as 0.04% in the figure caption and as 0.4% in the text.

I have access to Abaqus at present, but it is not equipped with a Fortran compiler. Thus, I am unable to modify and run simulations with subroutines.

I recommend:
a) delete the imposed DELTA_E as given in the subroutine V08.f, and replace it with the definition of eq. (6) in [1].
b) define straining of the model to 0.04% or 0.4% in V08.inp.
c) compare the output for the stimulus DELTA_E with Fig. 7.

If you do not succeed in reproducing the data from Fig. 7, then proceed to analyze V07.

As I was quite satisfied with my simulations listed unter 1) at this point and did not need the code any more, I discontinued that project at this point.

The zip-file contains:
a) all input-files. As the stimulus was defined in the subroutine, the geometry should be arbitrary, provided three directions in space exist as required by the definition of elasticity in the subroutine.
b) all Fortran-files. They differ only in their definition of the stimulus DELTA_E.
c) A Mathematica-Notebook 'Komarova-model.nb' for the analysis of 1) above, also printed to a pdf-file.

The attachment is a zip-file, but the site does not allow uploading zips. Change the file extension.

The above-mentioned Mathematica-Notebook also provides an analysis of the model in Mathematica, finding a numerical solution to the ordinary differential equations, in addition to the graphical representation of the Abaqus data for the simulations in 1) above. Both solutions generate a very good match to the data depicted in [1].

Abaqus was used in its release 6.9-EF1.

See also
http://imechanica.org/node/13153
for a list of publications that leak their codes in full or at least describe the algorithm in some detail.

You are free to use the code. Please kindly indicate the source of the code.

Kind regards

Frank Richter

References:
[3] and [4] were not studied while working on the project

[1] N. Bonfoh, M. Bilasse, P. Lipinski:
Modélisation du remodelage osseux
Revue de Mécanique Appliquée et Théorique, Vol. 1 (2008), no. 10, pages 717-726
 
[2] N. Bonfoh, E. Novinyo, P. Lipinski:
Modeling of bone adaptative behavior based on cells activities
Biomech Model Mechanobiol, vol. 10 (2011), pages 789–798

[3] S.V. Komarova, R.J. Smith, S.J. Dixon, S.M. Sims, L.M. Wahl:
Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling
Bone, vol. 33 (2003), pages 206–215

[4] S.V. Komarova:
Mathematical Model of Paracrine Interactions between Osteoclasts and Osteoblasts Predicts Anabolic Action of Parathyroid Hormone on Bone
Endocrinology, vol. 146 (2005), no. 8, pages 3589–3595