 |
|
ISSN 1320-0682 | ||||
| Volume 4 | January 1997 | ||||
Bruno Courtin, Anne-M. Perault-Staub, Jean-F. Staub
URA CNRS 1432, Faculte de Medecine Lariboisiere-Saint-Louis,
10 avenue de Verdun 75010 Paris, France
Email: courtin@ccr.jussieu.fr, staub@ccr.jussieu.fr
Numerical study of a reaction-diffusion model related to bone mineral metabolism reveals complex spatio-temporal behaviour in the vicinity of a Turing bifurcation. When these nonstationary patterns are explored in numerical experiments chosen for their biological relevance (size, geometry and growth of the two-dimensional reaction-diffusion space, random variations in parameter values), 3-D images and 2-D cross-sections of predicted spatial self-organisations can be produced as a time-record of the dynamic behaviour of the model with time being assimilated as an additional space dimension. Such images present striking analogies with the internal architecture of primary periosteal long bone. Our results suggest that dynamic properties of solute mineral species, intermediary between crystalline solid phase and solute parent ions, may play a major role not only, as previously reported, in the temporal self-organisation of bone mineral metabolism, but also in the emergence of the heterogeneous trabecular architecture typical of periosteal bone, at least during the first steps of embryonic bone development.
Most applications of nonlinear dynamic systems to morphogenesis and pattern formation in biology emphasise stationary "Turing" structures [37]. However, Meinhardt and Klingler [20] in their study of pattern formation on the shells of molluscs have used spatio-temporal behaviour such as travelling waves to explain the development of very complicated natural biological patterns. In the vicinity of a Turing bifurcation, nonstationary patterns resulting from a mixture of temporal instability and stationary spatial structure have been obtained numerically [27] or experimentally [18]. Bone mineral metabolism has been associated with a temporal [34] [35] [36] and spatial [4] [35] self-organised system. It has been suggested that some spatio-temporal behaviour of a reaction-diffusion model related to this metabolism could have qualitative analogies with the trabecular architecture of embryonic bone. Here, we go a step further with this model and specify numerical simulation conditions corresponding to a more realistic situation since it accounts, not only for geometry, size and growth of a diffusional space comparable to a given in vivo situation, but also for the existence of random fluctuations in parameter values.
Biological and physico-chemical
background
Bone development, as any process of embryonic tissue formation, is a continuum of sequential cellular decisions and expressions which are profoundly influenced by local extracellular and cellular interactions and which are responsible for giving this tissue its distinctive properties in terms of both morphology and physiological functions. While the bones in the body are of different shapes, lengths and embryonic lineages, one bone, the long bone of the embryonic chick (Gallus domesticus), has been chosen by biologists as a guide to document the common steps of osteogenesis. The first ossification is produced by the continuous supply of osteoblasts from the periosteum, a membrane which, in embryonic long bones, forms in the mid-diaphyseal region of the chondrogenic core. This periosteal ossification , intramembranous in nature, progresses by formation of a 3-D architecture of calcified trabeculae associated with vascularisation and stromal system development in the intertrabecular spaces. While its radial growth proceeds by an appositional process, periosteal bone expands proximally and distally as a sleeve over the growing cartilage core. Osteoblasts secrete an extracellular matrix (ECM) mainly composed of type I collagen plus a series of non-collagenous proteins (NCP) that further calcifies by precipitation of a Ca apatitic solid phase (HA) . The chick was chosen as an animal model for morphological analysis of periosteal long bone during its embryonic development because in this species, unlike mammals, there is no endochondral bone formation and thus no remodelling of metaphyseal bone until at least 7 days after hatching. As described by scanning electron microscopy [7], mineralisation of the chick tibia starts in ovo at approximately 8 days and, by 10 days, a cylinder of mineralised bone is clearly evident. At this time, scattered over its surface are regions of mineral deposition which are centres of growth. Later, these centres extend outward to form pillars that support a second mineral ring which has a diameter slightly greater than the first one. At 11 days, the periosteal bone sleeve is lengthened and the diaphyseal region is thickened by sequential deposition of additional cylindrical layers, the formation of one layer being not necessarily complete before the next begins. From 12 days to 2 days post-hatching, this process continues and the typical primary haversian structure (cortical bone) develops to form a concentric series of canals which radiate at an oblique angle with the axis of the diaphysis. This morphological description of periosteal bone development agrees with the pioneer work of Fell [9]. Using histochemical [28] and immunohistochemical [2] techniques, special attention has been more recently paid to the associated osteogenic cellular events and ingrowth of vasculature in intertrabecular spaces. This last process operates as follows: during the early events of osteogenesis, skeletal elements form in avascular positions; then, at the initiation of osteoid mineralisation, bone vasculature invades periosteum as precapillaries - that is, monolayered endothelial cells. At a further step (between 10 and 11 days of in ovo development), vascular invasion and erosion of the cartilage core occur; the nutrient artery penetrates the mineralised cylinder and the formation of bone marrow cavity begins, while the radial apposition of new periosteal bone proceeds continuously. Two days after hatching, a centrifugal circulation is established. Centrifugal movement of solute species in interstitial spaces of cortical bone [6] [15] [22] has been theoretically supported as largely dictated by both the difference in hydrostatic pressure between endosteum and periosteum and the porosity of the 3-D architecture of cortical bone with its particular vascular geometry [15]. However, the periosteum behaves as a peculiar region where solute nutrients accumulate near the periosteal surface, suggesting either some impermeability in this portion of bone or superimposition of the centrifugal endosteal circulation with the periosteal one, minimising in this way the influence of fluid flows on the transport of solute species in this bone region.
Under purely physico-chemical conditions, mineral precipitation in gel milieu is known for its ability to exhibit spatial [11] and spatio-temporal [39] patterns and this phenomenon was considered as a dissipative structure [25] associated with the dynamic properties of a reaction-diffusion system [14] [39]. We look here for the possibility that, during its development, the bone system uses such a 'generic' physical mechanism [24] by controlling both the reaction rates and the diffusivity of some mineral species involved in bone mineral metabolism. The mechanisms involved in the biological control of bone mineral homeostasis and matrix calcification are a matter of intensive research. In general terms, most current investigations focus on the capability of ECM constituents to regulate nucleation and crystal growth of bone solid phase(s) [1], [23]. Special attention was given to polyanionic NCP with their ability to form complexes with mineral, activating or inhibiting solid phase(s) production depending on their concentration and their state as matrix component, either immobilised or free in solution [19]. For instance, bone sialo-protein was proposed for initiating HA nucleation from ions at specific locations while osteocalcin is suggested to slow down HA crystal growth [30]. On the other hand, it is thought that, in vivo, diffusion has only local effects and that bulk diffusion plays no important role on the mineralisation process [3]. Indeed, although a concentration gradient of parent ions may exist in the immediate neighbourhood of a growing stable nucleus and may locally retard the growth or annihilate the appearance of other nuclei, both the homogeneity of Ca and phosphate ionic concentrations in extracellular fluids (the only solution-dependence considered in the current calcification schemas) and the limited thickness (some micrometers for osteoid) of the unmineralised extracellular spaces of bone under active osteoblasts are such that diffusion does not appear as a rate-limiting step.
Yet, a set of theoretical and practical considerations might lead to a substantial revision of this opinion. Firstly, in ionic solutions, not only free ions, but also an appreciable amount of semi-ordered and stable solute-solvent clusters are present [17]. While the classical scheme of nucleation and crystal growth postulates precrystalline aggregates as very unstable entities - that is, having a life span so short that their mutual interactions or their reaction with foreign species in heterogeneous milieu may be neglected - solute clusters have been proposed as actively participating in mineral precipitation [16]. Consistent with the existence of stable mineral species intermediary between free ions and the final solid phase in bone, a self-oscillating model of calcium metabolism exhibits the major role of two peculiar kinetic entities which are responsible for the dynamic circadian expression of bone mineral in vivo [34] and reminiscent of the existence in solution of low-ordered and semi-ordered solute clusters, respectively [26] [31]. Moreover, the various co-operative interactions between mineral clusters and bone matrix constituents must modify their stability, reactivity and effective and/or apparent diffusivity. Finally, diffusion in bone fluid occurs not only in the thickness, but also in the two other dimensions of the unmineralised bone matrix - that is, over a space whose size is a priori so large that, for particular reaction-diffusion rates, heterogeneity in the distribution of species like mineral solute clusters can be expected, associated or not with temporal self-organisation.
Our emphasis here is on the comparison between the sequence of the morphological events involved in embryonic bone formation and the features of spatio-temporal self-organisations generated by a reaction-diffusion model of bone mineral metabolism. An appropriate bridge between experimental and theoretical investigations could provide new perspectives on the co-operative mineral-organic-cellular interactions that contribute to the process of osteogenesis.
The model corresponds to the strictly nonlinear part of the temporal self-organised compartmental model of calcium metabolism in vivo [34] but, here, transport processes by diffusion were explicitly considered. Thus, the system of two-variable nonlinear partial differential equations is as follows:
where U and V denote the concentration of two species comparable to
solute clusters
mentioned above. U is formed from solute ions and V is a direct
precursor of bone solid phase. A is a constant flux of U formation from a
pool of parent ions maintained at a constant concentration. k1,
ka, kb, k2, k3
are non-negative rate constants defining respectively,
the dissociation of cluster U into parent ions, the spontaneous formation
of V from U, the nonlinear transformation (order-two autocatalysis) of U
to V, the dissociation of V into U and the use of precursor V for the
irreversible formation of solid phase.
DU and DV
are diffusion co-efficients for U and V
and 
, the Laplacian operator for each variable in a
two-dimensional (x,y) space:
with v = U or V, L the unit length, lx and ly being
the size of diffusional space
and the scale factors gx, gy in x and y dimension:
,
.
Geometry and growth of the diffusional space
When the first layer of differentiated osteoblasts secretes a ring of
matrix around the cartilage core of embryonic long bone, the bone
diffusional space is no thicker than the periosteum and may
approximate the surface of a cylinder if its thickness is neglected.
Consequently, boundary conditions were chosen as periodic and no-flux,
for x and y dimension, respectively. We used an initial cylinder of
radius r0=1.5 10-2cm and of length
ly(0)= 4.0 10-2cm.
The growth of this space was processed numerically by substituting lx for
and ly for
in Equation 2, where Gr
(6.94 10-8 cm.s-1) and Gl
(8.33 10-7
cm.s-1)
define the radial and the longitudinal (in both proximal and distal
directions of y dimension) linear growth rate, respectively.
Numerical resolution and simulation conditions
An explicit finite difference discretisation method with stepsize (dx,dy) not less than 0.5 10-3 and up to 1.0 10-3 cm was used, with automatic grid recalculation when a process of growth is explicitly taken into account. The ordinary differential equation system was integrated in double precision on a DECstation 5000/200, using the Gear method [12] with a relative error tolerance up to 1.10-12.
When the differential system was studied for its stable dynamic behaviour, the initial conditions were the randomly noisy (+/-1%) homogeneous steady-state value for U and V, and the numerical integrations were run over a time supposed to be long enough to reach the "asymptotic" solutions. When the behaviour in the presence of growth was studied, the spatial transient distribution of U and V was collected every simulation hour over the entire simulation duration. Two distinct situations were then considered:
The uniformly distributed random series is generated using the rand() function in C language.
Table 1 The set of used parameter values.
* specific values are given in legend figures.
Each figure was generated using Advanced Visualisation System
(AVS, Stardent Computer Inc). The 2-D figures were computed using
a blue- to red-coloured gradient between the minimal and median
(Vm) V values, the red colour representing
. Some
of them were obtained by computing longitudinal or transversal
cross-sections of the overall data used to generate 3-D
structures. 3-D images were depicted by computing the isosurface
of V between successive simulation times, with Vm as threshold value.
The calculated surface encloses the values of
. Under these
conditions, t is considered as a third dimension of space.
Cylindrical geometry of bone sleeve is obtained by transformation
of x dimension into polar co-ordinates.
The main results reported here are illustrated in Figure 1, Figure 2 and Figure 3.#Figure 3
Figure 1: 2D-Turing stationary structures self-generated by the model.
Figure 1 The 3 kinds of asymptotic stationary spatial behaviour for V concentration. 2-D Turing patterns have been generated for increasing R (
), (a) Hexagons (H1), R = 2.0 with DU=2.0 10-10, DV=1.0 10-10, (b) Stripes (S), R = 2.6 with DU=2.6 10 -10, DV=1.0 10-10, (c) Re-entrant Hexagons (H2), R = 6.0 with DU=6.0 10-10, DV=1.0 10-10. Spatial patterns are depicted as reported in Model and Numerical Methods; a red to blue gradient between median and minimal V value was used. A spatial stepsize dx=dy= 10-3 cm and a 90x40 grid was employed. Due to 'stable' defects, the structures do not appear completely regular.
Figure 2: Spatial self-organisation predicted by the model at a 2-D level comparable to cross-sections of the cylindrical bone sleeve.
Figure 2 2D-tranversal and longitudinal cross sections of 3D-structures such as that exemplified in Figure 3a. Slices were obtained by computing V values in (x,z) or (y,z) plane of the co-ordinate space used in AVS. The results have been obtained, for R=1.70, with DU=1.7 10 -10, DV=1.0 10 -10 and grid sized up to 220x280, in the presence (part b, d) or not (part a, c) of random fluctuations in parameter values. Mid-cylinder transversal sections (a) and (b) with cylinder radius of 3.5 10-2 cm and longitudinal sections (c) and (d) with a final cylinder length of 2.8 10-1 cm are shown. Spatial patterns are depicted as reported in Model and Numerical Methods; a red to blue gradient between median and minimal V value was used. Uniform grey parts illustrate the cartilaginous core.
Figure 3: Spatial self-organisation predicted by the model at a 3-D level comparable to the architectural arrangement of the cylindrical bone sleeve.
Part a: Overall view of the irregular external layer of the cylindrical sleeve.
Figure 3a 3D-image of spatio-temporal pattern obtained for R=1.70, with DU=1.7 10 -10, DV=1.0 10 -10 in the presence of random fluctuations in parameter values. The computed isosurface encloses the values of
. This image was obtained by using the AVS geometry viewer and the hardware rendering techniques (DECstation 5000/200). The overall view shows the irregular external surface observed after 44 h of simulated growth. The cylinder radius and length are 2.6 10-2 and 1.72 10-1 cm, respectively.
Part b: A detailed structural arrangement.
Figure 3b 3D-image of spatio-temporal pattern obtained for R=1.70, with DU=1.7 10 -10, DV=1.0 10 -10 in the presence of random fluctuations in parameter values. The computed isosurface encloses the values of
Part c: Another detailed structural arrangement.
Figure 3c 3D-image of spatio-temporal pattern obtained for R=1.70, with DU=1.7 10 -10, DV=1.0 10 -10 in the presence of random fluctuations in parameter values. The computed isosurface encloses the values of
. This image was obtained by using the AVS geometry viewer and the hardware rendering techniques (DECstation 5000/200). A detail of the internal architecture is given: arrangement between "pillars" and successive cylindrical layers shows the complex internal structure over a thickness of 2.0 10-2 cm in the mid-cylinder region.
In an earlier paper
[4], linear stability
analysis and numerical studies of the model (Equations 1 and 2) led us
to confirm the ability of this kind of nonlinear system to
generate spatial self-organisation. In spite of the temporal
instability of the homogeneous steady state associated with the set
of reaction co-efficient values (Table 1)
and responsible for the stable homogeneous solution corresponding
to a limit cycle (whatever the value of the diffusion
co-efficient ratio (
The assumption that
the geometry of the diffusional space
remains a cylinder and that both radial and longitudinal
growth is linear and isotropic, represents an idealisation of
the real system. However, it appears a reasonable simplification
since simulations are carried out only over a short duration
after initiation of the developmental process. Similarly, the
values of the initial space size and of the growth rates are
only approximate, although of correct magnitude if
compared to those of the biological system
[7] [28].
Random fluctuations in reaction co-efficient values have been
introduced explicitly during the simulations, mainly because
of the stability of the homogeneous limit cycle, eventually
co-existing with stable spatial or spatio-temporal behaviour for
sufficiently high R values. Indeed, as expected for a stable
behaviour, the choice of initial conditions belonging to the
stable limit cycle does not allow the reaching of the spatial and/or
spatio-temporal behaviour, even when initial conditions are
submitted to a very large noise. Although not shown here, a
detailed study of the effect of noise parameters on the dynamic
behaviour of our model indicates the existence, for
a reasonable (about +/-10%) amplitude, of a mean spatial and time
frequency of the random noise such that the limit cycle is destabilised
for the benefit of discernible, even if perturbed, stationary spatial
organisations.
Under such conditions, the range of R values for which spatio-
temporal organisation appears, increases particularly for the
low R values.
Figures 2 and 3 depict the kind of architecture predicted
by our model either at a 2-D level
(cross-sections of the cylindrical sleeve) or, as 3-D
images (overall and detailed views of external and internal
structure). However, since comparisons will be made with
histological
observations of a real bone system, we emphasise the two
following points to avoid any misinterpretation of
our results:
1) The assimilation of time, t, to one space dimension fits
well with the mode of appositional growth leading to the increase
in width of periosteal bone. It amounts to a time record of
the spatial distribution of one of the two variables.
Such a procedure that leads in this paper to 3-D
images, is quite analogous to that used by Meinhardt
[21] to establish similarities between
time records of a one-dimensional pattern and the 2-D spatial
organisations of mollusc shells. In our study, only the dynamic
behaviour related to the external surface of bone was considered.
Thus, the 3-D images generated from the model do not account for
internal bone remodeling, a process
responsible for modifications
of internal bone architecture through resorption and/or widening
of pre-existing calcified trabeculae. Consequently, direct
comparisons with real structures can be justified only under
conditions where remodelling can be neglected - for example, early after osteogenesis initiation.
2) The V variable is used to illustrate the spatial organisations
generated by our model. V is thought to be representative of the
diffusing semi-ordered solute clusters and considered, on the one
hand, as the mineral precursor for the HA solid phase precipitation
in extracellular bone matrix and, on the other hand, as reactive
entities involved directly or indirectly in cell-ECM
interaction (see last paragraph of this section).
Therefore, Figures 2 and 3 may
be interpreted only as the time record of the prepattern involved
in the embryonic bone-forming process.
Figure 2 illustrates the kind
of structural
organisation observed in two distinct situations,
according to the existence (Figures 2b and 2d) or not
(Figures 2a and 2c) of random fluctuations in parameter values
during the simulations. It visualises the transversal (in the
mid-cylinder region) and longitudinal cross-sections of the
cylinder for an R value of 1.70 - that is, lower than its critical
value, Rc - and after 80 hours of growth. At this time, the
thickness and length of the cylinder reach 0.2 mm
and 2.8 mm, respectively.
In both cases, spatio-temporal behaviour is obtained.
The temporal aspect of dynamics seems directly related to
the more or less continuous successive rings of high (red) and
low (blue) values of V concentration. The spatial aspect is
linked to the fact that, at some not regularly distributed
locations on these rings, high V values exist that ensure
connectivity between the successive rings. Note that, in the
absence of noise (Figures 2a and 2c), the size of low V value domains
enlarge when the radius and length increases, a result not observed in the
presence of noise. In this last case, because of starting from
initial conditions situated on the limit cycle, the earliest
behaviour is largely dominated by the temporal aspect, with only
few irregularities on the first ring. Later, the number of
these irregularities increases so that the structural complexity
is distributed over the whole space. Clearly, this kind of structural
organisation fits well with the histological observations of
periosteal long bone cross-sections, the high and low values
of V being associated with bone tissue trabeculae and
intertrabecular soft-tissue spaces, respectively.
Furthermore, the space between successive rings
(about 40 micrometers) and the size of the intertrabecular
spaces are not very different from the natural organisation
[7]
[28].
Figure 3a
shows an overall view of the external surface of the sleeve obtained
in the presence of noisy parameters and after 44 hours of growth - that is, at a time when the third ring is developing.
At a 3-D level, it appears as a complex surface which agrees
fairly well with the mineralised structure of embryonic chick
tibiae [7]. Irregularities in the thickening of
the cylinder result from the addition of an incomplete cylindrical
layer. Moreover, details presented in Figure 3b and
Figure 3c,
reveal amazing analogies between the predicted structure
and the process by which the
periosteal bone progresses under
physiological situations: from the second layer (Figure 3b),
singularities are born on a quasi-continuous smooth surface,
grow radially while increasing in the other dimensions and
then fuse to form a new 'continuous' layer. The spaces of low
V values - that is, the interconnected volumes situated between
the successive cylindrical layers supported by "pillars" (Figure 3c) - are regions available for
the localisation of the complex vascular network associated
with bone development. Although the model agrees well with the main events of the
developmental sequence of the first periosteal long bone
embryonic rudiments, obviously some details in the structural
organisation are not correctly reproduced: for instance, changes of
radial pillars supporting the successive cylindrical layers into
elongated "ridges" in the metaphyseal region together with the
further development of a regular concentric series of canals forming
the primary haversian structure
[7] are not
predicted here. However, for higher R values than those
used in this paper - that is, when H1 hexagonal stationary structure may be
obtained - there are similarities of the observed pattern with that
simulated, including changes in
the angle that the haversian canals form with endosteal surface
[6].
This could suggest an overall patterning development as follows:
At the initiation of osteogenesis, the reaction-diffusion domain
is so small that the temporal organisation predominates with the
formation of a first homogeneous cylinder and only few irregularities
scattered over its surface; when a sufficiently large diffusional
space is reached, the spatial organisation bifurcates towards
a stationary hexagonal pattern generating the typical haversian
structure. Between these two stages, there is transient behaviour
with prevalence of spatio-temporal dynamics like that illustrated
in this paper.
Note that, according to the model, such a sequence from homogeneous
temporal to stationary spatial behaviour is obtained only if the
system is submitted to sufficiently high random fluctuations
in parameter values.
Being aware of Wolpert's expert assertion
[38], 'do not infer process from pattern,
since so many processes can produce the same pattern', we do
not expect to justify the series of assumptions and
simplifications included in our model, on the sole "satisfactory"
comparison between predicted patterns and observed bone spatial
structures. However, at least one point that we regard as thoroughly
important seems on the way to being reached: a single dynamical
process might have a major role in both temporal and spatial
self-organisation of bone mineral metabolism and
this connects two aspects usually studied separately:
firstly, the mineral homeostatic function of bone
[23] [29]
and secondly, the development of an embryonic architecture involved
in the further structural adaptation
of bone to mechanical usage
[5] [10].
Yet, at the very least, some important questions about the
interpretation of our model have to be answered. They mainly
rest on biological considerations and shall be largely argued
elsewhere [33]. Briefly, they concern:
We thank L.Toubiana for his decisive contribution in the first stages
of this study and V. Myrtil for his technical assistance.
We are indebted to the Centre de Calcul Recherche et Reseau Jussieu (CCR Paris 6 et 7)
for the computation time allocated. We are also
grateful to Ph. Depondt for improvement in the English style of the manuscript.
This work was supported by Centre National de la Recherche Scientifique (CNRS).
Stationary and non-stationary patterns
in a 2 dimensional R-D space
)),
a Turing bifurcation - that is,
potential stationary spatial structures - occurs for a critical
value of R (
). As illustrated in
Figure 1, the various well-known
stationary 2-D spatial structures
[8], hexagons-type 1 (H1), stripes (S) and
re-entrant hexagons (H2), can be obtained for increasing R values
beyond Rc and for a size of diffusional space corresponding
to the initial surface of the cylinder. Note that, contrary to the
usual situation, there is no requirement for a high value of
Rc - that is, for very different diffusion co-efficients
[27]. Moreover, nonstationary spatial
patterns associated with spatio-temporal chaotic dynamics
have been shown to exist in the immediate vicinity
of Rc, due to the interaction of the Turing bifurcation with the
temporal instability. In this paper, we have
chosen to study the type of behaviour that presents the advantage of
not being mutually exclusive of temporal and/or spatial
self-organisation.
The basic hypothesis is that this behaviour might be applied
to periosteal long bone development in a chick embryo.
Thus, we have
attempted to take into account both the geometry and growth
of this biological system as well as the fact that, like any physical system, this
biological system is submitted to stochastic perturbations.Geometry, growth and random fluctuations in
parameter values
Time recording of non-stationary behaviour on a cylinder
and trabecular bone architecture
3-D overall and detailed views
References
Glossary
Diagrammatic representation of a post-embryonic long bone.
Complexity International (1997) 4