UNIVERSIDAD COMPLUTENSE DE MADRID 
 FACULTAD DE CIENCIAS MATEMÁTICAS 
DEPARTAMENTO DE MATEMÁTICA APLICADA 

 

 

  
 
 

TESIS DOCTORAL      
 

Modelos matemáticos de reparación ósea 
Mathematical models of bone repair 

 
 

MEMORIA PARA OPTAR AL GRADO DE DOCTOR 
 

PRESENTADA POR 
 

    Luis Echeverri Delgado  
 
 
 

Directores 
 

Miguel Ángel Herrero 
Gerardo Oleaga 

 
 
 
 

Madrid, 2014 
 
 
 
 
 
©Luis Echeverri Delgado, 2014 



MODELOS MATEMÁTICOS DE REPARACIÓN ÓSEA
MATHEMATICAL MODELS OF BONE REPAIR

Memoria presentada para optar al t́ıtulo de Doctor en Ciencias
Matemáticas

Por:

Luis Echeverri

Dirigida por:

Dr. Miguel Ángel Herrero
y

Dr. Gerardo Oleaga

Facultad de Ciencias Matemáticas

Departamento de Matemática Aplicada

UNIVERSIDAD COMPLUTENSE DE MADRID

2014



ii

Universidad Complutense de Madrid



Luis:-¿y tu sabes a quien quiero yo más?-
Ana (3 años de edad):-“a Yo, a Yeye, y a la Abuela.”-

iii



iv

Universidad Complutense de Madrid



Contents

Resumen 1

Abstract 5

1 Preliminaries 9
1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2 Bone structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3 Bone remodelling and the Basic Multicellular Units . . . . . . . . . . 13
1.4 Fracture Healing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.5 Mathematical models . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

1.5.1 Mechano-regulatory models . . . . . . . . . . . . . . . . . . . 19
1.5.2 Bio-regulatory models . . . . . . . . . . . . . . . . . . . . . . 20

Bibliography 23

2 Bone remodelling 29
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.2 A minimal software for BMU operation . . . . . . . . . . . . . . . . 32

2.2.1 The role of molecular signals in bone remodeling . . . . . . . 32
2.2.2 Signal integration by individual cells. Inhibitory proteins . . 34

2.3 Mathematical modelling . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.3.1 Cell algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.2 Spatial dependence in the model . . . . . . . . . . . . . . . . 41
2.3.3 Cell movements . . . . . . . . . . . . . . . . . . . . . . . . . . 41

2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.4.1 The BMU software proposed successfully recapitulates bone

remodeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
2.4.2 The size of the BRC adapts to the number of apoptotic OCY 43
2.4.3 Osteoclast lifespan depends on the depth of the BRC . . . . . 44
2.4.4 The BMU software is robust to signal noise . . . . . . . . . . 44

v



CONTENTS vi

2.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

Bibliography 49

3 Early stages of Bone fracture healing 53
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.2 Basic facts on bone fracture healing. The first 48 hours . . . . . . . 54
3.3 A mathematical model for the formation of a fibrin-collagen callus

template. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.3.1 Generation of fibrin monomers upon thrombin-induced fibrino-

gen cleavage. . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.3.2 Estimating fibrin clot formation . . . . . . . . . . . . . . . . . 61
3.3.3 Platelet activation and growth factor release . . . . . . . . . 64
3.3.4 Recruitment of mesenchymal cells (MSC). Formation of early

fibrin-collagen callus (EFC) . . . . . . . . . . . . . . . . . . 65
3.4 A sequence of times . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
3.5 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Bibliography 75

A Details of computer simulations of the model of BMU operation
during bone remodeling 81
A.1 Time evolution of the BMU . . . . . . . . . . . . . . . . . . . . . . . 81

A.1.1 Pseudo-code guidelines . . . . . . . . . . . . . . . . . . . . . . 85

B Experimental procedure 89

Universidad Complutense de Madrid



Resumen

En esta memoria se obtienen y analizan algunos modelos matemáticos de remod-
elación y reparación ósea. Para ello, y tras un primer caṕıtulo introductorio en el
que se presentan resultados preliminares para los estudios posteriores, se aborda en
el Caṕıtulo 2 la modelización del mecanismo de mantenimiento que tiene lugar a
lo largo de la vida de cada persona, y en virtud del cual en pequeñas regiones del
esqueleto el hueso viejo es reemplazado por el nuevo de manera que la cantidad total
de hueso permanece constante.

En concreto, en el caṕıtulo 2 se formula un modelo mı́nimo para el funcionamiento
de las llamadas BMU (Basic Multicellular Units), grupos de células diferenciadas
que actúan coordinadamente para eliminar hueso en regiones marcadas para ello
(generalmente mediante señales emitidas por osteocitos, que actúan como sensores
mecánicos) y a continuación generan una matriz extracelular que, tras su mineral-
ización, da lugar al nuevo hueso que reemplaza al antiguo. La actuación de tales
unidades orgánicas es conocida desde los años 60 del siglo pasado, y aunque su
esquema general de actuación es conocido, y muchas de las señales qúımicas in-
volucradas han sido identificadas, el mecanismo preciso por el que tales unidades se
forman cuando es necesario y se deshacen una vez cumplida su tarea, sigue sin ser
bien conocido.

En esta memoria proponemos un modelo operativo simple para BMUs, que se
basa en las siguientes hipótesis. En primer lugar, cada célula en una BMU puede
elegir una entre un número muy limitado de posibles acciones a realizar (dividirse,
diferenciarse a otro tipo celular o morir). Tales decisiones no son aleatorias, sino
que estn reguladas por inhibidores internos espećıficos que bloquean cada una de
esas acciones, de modo que el primero de ellos cuya concentración baja de un cierto
valor cŕıtico determina de forma irreversible la acción a seguir. La dinámica interna
de estos inhibidores está modulada por señales externas, provenientes en general de
otros tipos celulares presentes en la BMU. Mostramos en este caṕıtulo que basta con
un número muy pequeño de tales señales (tres) para que el proceso resulte operativo.
En cualquier caso, cada célula decide individualmente lo que hace en cada momento,
sin que haya ninguna determinación previa para ello, ni genética ni de otro tipo.

1



CONTENTS 2

Sin embargo, el resultado de esta suma de decisiones individuales independientes
es una acción coordinada y perfectamente regulada, que permite remodelar la zona
elegida para ello de manera robusta y eficaz, como se observa mediante el análisis
de los resultados obtenidos al simular el modelo correspondiente, y que se presentan
en el mismo caṕıtulo. Los resultados obtenidos permiten formular cuestiones rela-
cionadas con la presencia detectada de un número de señales mucho mayor que las
estrictamente necesarias, lo que probablemente está relacionado con conceptos como
la multifuncionalidad (una misma señal induce diversos efectos en tipos celulares
distintos) y la redundancia (en el sistema considerado persisten diversos mecanismos
para alcanzar un resultado, que garantizan el funcionamiento del sistema cuando
uno falla, pero que representan un coste metabólico para las células involucradas).

El esquema que acabamos de resumir ilustra un aspecto de la metodoloǵıa seguida
en esta memoria, que también se usa en el siguiente caṕıtulo, y que nos parece par-
ticularmente relevante. En concreto, nuestro objetivo en cada uno de los problemas
estudiados es obtener modelos tan simples como sea posible, que permitan reproducir
los hechos observados y más aún, que proporcionen indicios sobre el funcionamiento
interno de las estructuras consideradas que no han sido previamente incluidos en la
formulación del modelo. Este enfoque es distinto de la estrategia habitual de mod-
elización, que busca incluir tantos datos biológicos como sea posible en cada modelo,
lo que conduce a una gran complejidad en los sistemas de ecuaciones obtenidos y
en particular hace precisa la introducción de gran número de parámetros, sobre los
que se sabe muy poco. Ello obliga a elegir tales parámetros para ajustar cada ex-
perimento concreto que se analiza, sin que se obtenga aśı un conocimiento real de la
dinámica involucrada, ni se pueda explicar de manera convincente las razones de los
comportamientos observados.

En el tercer caṕıtulo de esta memoria se estudia la sucesión de procesos que tienen
lugar en las etapas tempranas de la reparación del hueso después de producirse una
fractura, y que son cruciales para una completa curación del mismo. Se proponen
modelos matemáticos que permiten estimar los tiempos caracteŕısticos de tales pro-
cesos en términos de sus parámetros cinéticos. Para ello, se analizan tales etapas
de manera separada, hasta la formación del primer callo fibroso, que proprociona el
molde básico sobre el que posteriores procesos de remodelación (no estudiados en
esta memoria) darán finalmente lugar a nuevo tejido óseo, plenamente funcional.

El modelo propuesto es de tipo modular, de modo que consiste en una con-
catenación de modelos simples que se suceden uno a uno hasta completar la etapa
temprana de reparación estudiada. Se comienza por describir la formación de un
coágulo sangúıneo inducido por la fractura, cuya aparición se caracteriza como una
transición de fase tipo sol-gel, en la que el coágulo resultante es descrito como un
gel. La formación de dicho coágulo, que debe conectar los bordes de la fractura
proporcionando una primera unión entre ambos, es condición necesaria para el éxito
del proceso completo de reparación . Describimos en primer lugar la generación
de monómeros de fibrina a partir de la activación de su precursor inactivo presente
en la sangre (fibrinógeno) mediante la acción de trombina, liberada tras las roturas
vasculares inducidas por una fractura. Los monómeros aśı formados se polimerizan

Universidad Complutense de Madrid



CONTENTS 3

rápidamente para formar grandes agregados. Usando la teoŕıa de polimerización
clásica, el proceso correspondiente puede ser modelado por un sistema infinito de
ecuaciones de coagulación-fragmentación para los k-meros con k ≥ 1, en las que
los coeficientes de coagulación son proporcionales al número de enlaces funcionales
libres en cada monómero. Dado que no es posible en general resolver las ecuaciones
anteriores, se recurre a la introducción de una relación de clausura para caracterizar
la transición de fase buscada. Ello se hace sobre las ecuaciones para los momentos de
la distribución de poĺımeros, en las que ahora es posible estimar de manera simple
el cociente M2/M1, que corresponde al tamaño medio molecular, de modo que la
transición de fase buscada aparece cuando este cociente desarrolla una singularidad.
Se obtienen fórmulas para el tiempo de coagulación que son expĺıcitas en varios casos
y que permiten estimar el tiempo de coagulación en el caso más general.

Por otra parte, además de romper las moléculas de fibrinógeno, la trombina activa
a las plaquetas de la sangre, que cambian su forma liberando factores de crecimiento.
En este caṕıtulo se analiza el efecto inducido por un factor de crecimiento genérico
g, que se difunde a través del coágulo para atraer a las células mesenquimales. Éstas
juegan un importante papel en la formación del primer callo elástico, analizado en
este trabajo, y que proporciona el primer molde con estructura elástica suficiente
para servir de base a posteriores procesos de osificación no considerados en esta
memoria.

Tras proponer un modelo de tipo clásico para la difusión de g, se obtiene una
solución expĺıcita para dicha función, cuyo gradiente provoca el movimiento por
quimiotaxis de las células mesenquimales m desde el borde la fractura hacia su
interior. Se muestra que el perfil del gradiente de g se estabiliza rápidamente y
puede aproximarse por una función lineal, siempre que el intervalo temporal no
supere al tiempo de vida media de las plaquetas. Se obtiene aśı una estimación para
el tiempo TF de formación del callo fibroso, que completa las obtenidas previamente
para otros tiempos caracteŕısticos del proceso. Los resultados aśı obtenidos son
finalmente comparados con los existentes en la literatura.

Universidad Complutense de Madrid



CONTENTS 4

Universidad Complutense de Madrid



Abstract

This work is concerned with the mathematical modelling and simulation of bone
remodelling and repair. To this end, after a first introductory chapter where a
number of preliminaries are gathered, we consider in Chapter 2 the modelling of
a homeostatic maintenance mechanism that is at work during a whole life, and by
means of which in small, targeted regions in the skeleton old bone is first destroyed
and then replaced by new one, so that the total amount of bone remains constant.

More precisely, in Chapter 2 a model for the operation of the so-called Basic
Multicellular Units (BMU) is formulated. BMUs are groups of differentiated cells
that act co-ordinately to destroy (resorb) bone in targeted zones, identified as such by
osteocytes, that are able to act as mechanical sensors, and then secrete extracellular
matrix which upon subsequent mineralization yields new bone to substitute the
previous one. The existence of such organic units is known since the sixties of
the former century and, while its general way of action is known, and many of the
chemical messengers involved have been identified, the precise manner by which such
units are formed when needed, and disbanded afterwards, remains to be ascertained
as yet .

In this memoir, a simple operative model for BMU operation is proposed. It is
based in the following assumptions. To begin with each cell in a BMU can only select
one among a very limited choice of actions (to divide, to die or to differentiate into
another cell type). Such decisions are not randomly made, but are instead regulated
by specific internal inhibitors, which block each of these actions, so that the first of
them that goes below a critical threshold irreversibly determines cell fate. The inter-
nal dynamics of such inhibitors is modulated by external signals, released by other
cell types present in a BMU. We show that a very reduced number of such signals
(actually three) suffices to keep BMU at work. In any case, each cell individually
decides its course of action without resorting to any previously determined plan, ge-
netic or otherwise. Nevertheless, the outcome of this sum of independent individual
decisions turns out to be a coordinated and perfectly regulated course of action, that
leads to a robust and efficient remodelling of the zone involved, as can be seen from
analysis of the results obtained upon simulation of the model, that bare described

5



CONTENTS 6

in the same chapter. The results thus obtained lead to interesting questions related
to the detected presence of a set of signals much larger than that actually needed
for the model. This is probably related with aspects as multi-functionality (a same
signal may induce different actions in different cells) and redundancy (coexistence in
a same system of different mechanisms that yield the same result, thus ensuring a
security mechanisms against failure, though at the expense of a metabolic costs for
the cells involved).

The pattern described above shows one of the methodology characteristics fol-
lowed in this thesis, which is also used in the next chapter, and that we consider
especially relevant. Namely, our goal for each of the problems under study is to
obtain simple models which can recreate the observed facts. Furthermore, these
simple models will provide indications on the internal behaviour of the considered
structures, which had not been included previously in the formulation of the model.
This approach is slightly different to the usual modelling strategy that is to include
as much detail as possible, leading to typically very complex systems of equations
with a large number of parameters, most of which are unknown. This forces one
to identify such parameters for each specific analyzed experiment, without gaining
any real knowledge on the dynamics of the problem or explanation of the observed
behaviour.

In the third chapter of this thesis the processes which take place in the early
stages of a fractured bone repair are studied. These early stages are crucial for the
complete recovery of the bone. Mathematical models are provided to estimate the
time scales of these processes as a function of their kinetic parameters. To do so
the different stages are analysed separately, up to the formation of the first fibrous
callus, which is the basic mould on which further remodellation processes will take
place (not studied in this work). This process ends with fully functional new bone
tissue.

The proposed model has a modular structure, formed by a sequence of simple
models which follow one another until completing the early stages of the bone repair
under study. Initially, the onset of a blood clot induced by the fracture is described as
a phase transition of sol-gel type, where the clot is represented as a gel. The formation
of this blood clot, which connects the edges of the fracture, is a necessary condition
to ensure the repair process is successful. The first stage is the generation of fibrin
monomers from the activation of their inactive precursor present in blood (fibrinogen)
by the action of thrombin, which is released after the vascular breaks induced by a
fracture. The monomers which have just been created polymerise rapidly to form
large agglomerates. Following classical polymerisation theory, the above process can
be modelled using an infinite system of coagulation-fragmentation equations for the
k-mers with k ≥ 1. The coagulation coefficients are proportional to the number of
free functional bonds of each monomer. In general the mentioned equations cannot be
solved, thus the introduction of a closure relation to characterise the phase transition
is needed. This is done for the equations describing the moments of the distribution
of polymers, for which the coefficientM2/M1 can be easily estimated. This coefficient
describes the molecular mean size, therefore the phase transition will be present when

Universidad Complutense de Madrid



CONTENTS 7

a singularity appears in the coefficient. Explicit formulas for the coagulation time
are obtained for several cases, and these provide an estimate for coagulation time
for a more general case.

Thrombin, apart from breaking fibrinogen molecules, also activates blood platelets,
which change shape whilst releasing growth factors. In this chapter the effect induced
by a generic g growth factor is analysed, which spreads around the clot to attract
mesenchymal cells. These play an important role in the formation of the first elastic
callus (analysed in this work) that provides the first elastically structured mould.
This mould is the building foundation for posterior ossification processes which are
not considered in this work.

A classical diffusion model for g is solved, obtaining an explicit solution. The
gradient of g induces a chemotaxis movement of the mesenchymal cells m from the
edge of the fracture towards the interior. The profile of the gradient of g is shown to
stabilise rapidly and can be approximated by a linear function, as long as the time
interval does not overcome the half-life of the platelets. In this way an estimate of
TF is obtained, which is the time of formation of the fibrous callus. This completes
the estimates obtained previously for other characteristic time scales. The obtained
results are compared to those found in the literature.

Universidad Complutense de Madrid



CONTENTS 8

Universidad Complutense de Madrid



Chapter 1
Preliminaries

1.1 Introduction

The human skeleton is a complex structure which provides a supporting framework
for the whole body. In addition, it acts as a shield to protect internal organs, serves as
a reservoir for a number of ions (particularly for calcium) and plays a crucial role in
body motion by mediating the transmission of force arising from muscle contraction.
Bone function is thus essential for human life, which makes important to keep bone
homeostasis by preventing, and when needed repairing , all kinds of bone disorders.

Out of the many bone-related hazards, bone fracture is a particularly relevant
one. About 1.5 million people suffer a bone disease-related fracture each year in the
U.S [54], and recent studies claim that about 280 fragility-related fractures per hour
occurred in Europe’s five largest countries and Sweden in 2010 [52]. Such figures are
expected to keep increasing in the future, particularly in aging societies as developed
countries are.

In spite of such sobering data, it turns out that human bone has an impressive
capability for self-repair. In fact, our skeleton is continuously being renewed, so
that old bone is being destroyed in some places, to be eventually replaced by new
bone produced in the same locations. In that manner the whole skeleton is changed
every 10 to 15 years [37]. This self-renewal process, termed as bone remodeling, is
known to be triggered by local changes in mechanical load, and takes place at a
microscopic scale. Such repairing task is carried out in each case by a comparatively
small number of cells, that make up a Basic Multicellular Unit (BMU) [18] which is
mobilized for that purpose and then disbanded once its task is concluded. In fact,
BMUs are recruited when osteocytes, a characteristic bone cell type, detect that
mechanical stress goes locally outside of a range of admissible values [43].

A particular case in which the previous situation applies corresponds to the onset
of micro-fractures . These can be induced just by physical exercise (including walk-
ing ) and their occurrence (and efficient repair) is believed to significantly contribute

9



CHAPTER 1. PRELIMINARIES 10

to keep bone remodeling in good shape [49]. However, when many micro-fracture
coalesce, or when a comparatively large fracture appears (for instance , as a conse-
quence of a traumatic event), BMUs turn out to be too small to take care of the
repairing process then needed , which has to be carried out at larger space and time
scales. To address the resulting new challenges, a number of strategies have been set
in place in the course of evolution. In particular, inflammatory signals play now a
key role and a blood clot is formed that , upon serving as a first bridge to connect
fracture borders, is subsequently remodeled to produce several callus with increased
stiffness, which will be eventually replaced by new bone when the process successfully
concludes. It is worth to be noticed that, unlike other tissue healing processes that
leave scars when concluded, bone remodeling and bone fracture repair are scare-less,
so that no interface between old and new bone can be easily tracked. This seems
to be related to the fact that bone healing shares several characteristics in common
with the developmental program according to which ossification is carried out at the
embryonic stage [20].

In this memoir we shall address two problems, each related with one of the two
repairing events described before. To begin with, in our next Chapter 2 we will
provide a simple mathematical model that accounts for successful BMU operation
during standard, homeostatic bone remodeling. In that Chapter, BMU operation
will be described as an emergent behavior that can be derived from a reduced set
of assumptions on the possible individual choices of the cells involved in one hand,
and on the number of effective chemical cues mediating such decisions on the other.
Then in Chapter 3 we shall be concerned with fracture repair, in particular in the
case of fractures small enough so that theoretical results could be compared with
experimental ones, but beyond the range of BMU operation alone. In particular,
the onset of a blood clot will now be crucial and has to be taken into account in the
corresponding model. It will be shown there that preliminary (albeit in our opinion,
illuminating) estimates on the characteristic times of the repairing mechanisms can
be deduced from a simple set of assumptions concerning kinetic aspects of the early
stages of the corresponding bone healing process.

We have used several times the term simplicity when describing the goals ad-
dressed in this work. Indeed, a basic pattern in our approach consists in showing
that key aspects of the processes considered can be retrieved from the dynamics of
mathematical models where the underlying molecular and cell mechanisms have been
reduced to their minimal possible ingredients. No attempt to build comprehensive,
exceedingly complex models requiring of extensive parameter fitting is made here.
We have rather chosen a different path, and in each case have tried to ascertain the
minimal number of mechanisms required to account for the processes considered.

The results previously sketched (in a very succinct manner) are to be found in
Chapters 2 and 3, as already mentioned. The rest of this introduction is devoted to
provide a number of preliminaries which will be later referred to in such Chapters.
In particular, we will shortly recall a few facts about bone structure in Section 1. 2
below. Additional information about bone remodeling and the role therein played
by BMUs can be found in Section 1.3, whereas a brief review of facts concerning

Universidad Complutense de Madrid



1.2. BONE STRUCTURE 11

bone fracture will be gathered in section 1.4. The emphasis in these sections is
on the biological facts underlying such processes. A short discussion on the use of
mathematical models to better understand the mechanisms involved is then provided
in Section 1.5.

1.2 Bone structure

The skeleton in vertebrates is formed by two tissues: bone and cartilage. In humans
the skeleton contains over 200 different elements that have been traditionally clas-
sified according to shape into long, short, flat and irregular bones. Cartilage is the
conective tissue that can be found in the joins between bones. It is formed by cells
called chondrocytes that produce a large amount of extracellular matrix of collagen
fibers. One of the main differences between bone and cartilage is that bone is a
highly vascularized tissue and cartilage has virtually no blood vessels [51]. In this
work we are going to focus on bones and we will not study cartilage in deep.

The individual bones could be described at the macroscopic and/or microscopic
scales. At the macroscopic scale cortical and trabecular bone can be distinguished.
The cortical or compact bone is located in the outer surface of the bones, is harder,
stronger and stiffer than cancellous bone and contributes to the 80% of the weigth of
human skeleton [37]. The trabecullar or cancellous bone is found in the inner part of
flat bones, is a highly porous, open-celled structure formed by interconected plates
and needle-like structures (named trabecullae) of bone tissue. The intertrabecullar
space is filled with marrow containing bloood vessels, nerves and several cell types.
In long bones the cancellous bone is located at the ends while in short bones or flat
bones are located everywhere [37].

Figure 1.1: Bone macro-structure of the upper extremity of left femur. (Re-
produced according to the Creative Commons Atribution 3.0 Unported Licence
http://creativecommons.org/licenses/by/3.0/legalcode)

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 12

At a microscopic scale bone has two different structures, lamellar and woven bone.
Lamellar bone is formed by concentric parallel layers of crystals and mineralized
hardened collagen in a system of short structures known as osteons or Haversian
system. This bone is a highly organized, anisothropic, and slowly formed material.
The woven bone is a fast formed and poorly organized material in which collagen
fibers and crystals are more or less randomly distributed [37]. Most of bone tissue in
adult is lamellar, woven bone appearing mainly as a consequence of fracture healing.

Figure 1.2: Bone micro structure. Osteons are roughly cilindrical structures that are several
millimeters long and around 0.2 mm in diameter [17]. (Reproduced according to the Creative
Commons Atribution 3.0 Unported Licence http://creativecommons.org/licenses/by/3.0/legalcode)

Several cell types are present in the main stages of bone formation: osteo-
clasts, osteoblasts, bone lining cells and osteocytes. Bone lining cells, also called
pre-osteoblasts or surface osteoblasts, can be found in the periosteal and endosteal
surfaces of bones, whereas osteocytes (differentiated osteoblasts) remain trapped in
the interior of the mineralized bone forming cavities named lacunae. Osteoblasts,
or bone forming cells, are derived form mesenchymal stem cells (as are bone lining
cells and osteocytes) and are located in the periosteum. The main product released
by osteablasts is type I collagen, but osteoblasts also secrete other noncollagenous
proteins. Osteoblasts form bone by facilitating mineralization but the mechanisms
by which this occur are not well understood [14].

Osteoclasts, or bone resorbing cells, are derived from hematopoietic stem cells.
Its main function is to destroy (resorb) mineralized bone and calcified cartilage. An
osteoclast is a multinucleated cell, and human osteoclasts on bone have five nuclei
and are about 150-200 µm in diameter [22]. (see Figure (1.3))

Bone develops by the process of ossification or osteogenesis. During this process
osteoblasts secrete osteoid, an amorphous material that gradually becomes densely
fibrous. Calcium phosphate crystals are deposited in the osteoid thereby becoming

Universidad Complutense de Madrid



1.3. BONE REMODELLING AND THE BASIC MULTICELLULAR UNITS 13

Figure 1.3: Bone cells that are present in the main stages of bone forma-
tion. (Reproduced according to the Creative Commons Atribution 3.0 Unported Licence
http://creativecommons.org/licenses/by/3.0/legalcode)

bone matrix. This last part of the process is termed mineralization. Two different
modes of osteogenesis are observed: intramembranous ossification, in which bone
formation occurs directly in primitive highly vascular connective tissue, and endo-
chondral ossification, in which bone formation takes place in pre-existing cartilage.
Endochondral ossification occurs in long bone formation (see Figure (1.4)) and in
fracture and osteotomy repair [13]. The terms intramembranous and endochondral
refers to the type of tissue being replaced, not to the type of bone synthetized, which
is the same in both cases [50].

Bone growth and remodelling are two different processes in the skeleton of verte-
brates. In the first one bone is formed in sites where it was not present before, results
in changes in shape and/or architecture of the bones, and occurs since the first stages
of fetal development until the second decade of life in humans [37]. The second pro-
cess, i.e. replacement of old bone by new one and reparation of microdamage [43],
takes place throughout the whole of life and, as we have already remembered, is re-
quired for mineral homeostasis, for adapting to mechanical changes and to maintain
the healthy functionality of bones.

1.3 Bone remodelling and the Basic Multicellular
Units

We have mentioned that bone growth and remodelling are two different processes, a
main difference being the destruction of previously existing bone. The replacement
of old bone by new one is a process that involves resorption and formation of bone

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 14

Figure 1.4: Formation and growth of long bones by endochondral ossification. Mesenchymal
condensation leads to a cartilage ”model”. Chondrocyte undergo apoptotic and cartilage mineral-
ization starts. Vascular invasion from vessels allows the migration of osteoblast precursor cells that
deposit bone on the degradated matrix scaffold. (Reproduced according to the Creative Commons
Atribution 3.0 Unported Licence. http://creativecommons.org/licenses/by/3.0/legalcode)

without participation of a third tissue (i.e. case of endochondral ossification). Re-
modelling required for mineral homeostasis need not to occur in specific sites and
remodelling related to mechanical changes and damage is site-specific.

The resorption of bone is carried out by osteoclasts acting on small pockets
of bone and followed by osteoblasts precursors recruitment that differentiate and
replace the resorbed bone by a suitable extracellular matrix (ECM) where bone
will be eventually produced. Resorption and formation processes are coordinately
carried out by discrete temporaly anatomic units called ”Basic Multicellular Units”
(BMUs) first described by Harold Frost in the sixties [21]. The region where BMU
operates, referred to as the bone remodelling compartment (BRC), is a canopy of
bone lining cells/osteoblasts and a nearby capillary covering the BMU [49].

BMU are present in cortical and trabecular bone and differ in structure and in
the ways they remove and replace bone. In trabecullar bone, BMU are located at
the surface covered by a canopy of cells of mesenchymal origin. The resorbing zone
is cleaned by lining cells and macrophages and osteoblasts precursors differentiate
to fill the gap created by osteoclasts. BMU in the cortex display a moving team of
osteoclasts in it’s front forming the cutting cone or hemicone, the remaining space
being filled by blood vessels, nerves and conecting tissue. (see Figure 1.5).

BMU move in three dimensions excavating and refilling a tunnel in cortical bone
or a trench across the surface of cancellous bone. A cortical BMU travels about 4000
µm at about 20 µm/day, taking about 200 days to complete it’s task. A cancellous
BMU travels about half this distance about half the speed, taking about the same
period of time to fulfill that process [43]. The resorption period has a duration

Universidad Complutense de Madrid



1.3. BONE REMODELLING AND THE BASIC MULTICELLULAR UNITS 15

Trabecular bone

Cortical bone

cement line

New bone

old bone

osteoid

blood vessel

osteoblasts

osteoclast

osteoblasts

osteocytes

osteoclast

Intracortical surface 

of Haversian canals

trabecular 

surface

precursor cells

Macrophage

Osteoclast precursor Osteoblast

precursor

T cell

Osteoclast 

precursor

Figure 1.5: Remodeling is initiated within BRCs at points beneath the canopy of cells lining
trabecular bone (upper panels) and within cortical bone Haversian canals (lower panels). Osteo-
clasts (OCs) are formed from hemopoietic precursors (HSC) supplied by marrow and the blood-
stream. Precursors of osteoblasts come from MSCs in the marrow, from blood and from pericytes,
and differentiate within the BMU through the osteoblast precursor stage to fully functional syn-
thesizing osteoblasts and to osteocytes; lining cells may also differentiate into active osteoblasts.
T cells and macrophages can gain access to the BRC from the blood supply. (Reprinted by
permission from Macmillan Publishers Ltd: BoneKey Reports 3. Article number: 481 (2014)
doi:10.1038/bonekey.2013.215, copyright 2014)

between 30-40 days and is followed by the formation period with a mean duration of
150 days [18]. Remodelled regions can be recognized by cement lines that follow the
festoon-shaped surfaces that were eroded by osteoclasts [49]. In normal bone, the
result of this process is the complete refilling of the resorption lacuna. The succesful
coordination of bone resorption and formation in each BMU is refered to as coupling
(see Figure 1.6). The mechanisms responsible for such coupling remain largely elusive
[18]. Any disturbance in coupling may result in pathologies as Osteoporosis, which
appears when osteoblasts are not able to fully refill the resorption region, so that a
net loss of bone occurs in each remodelling event. Conversely, Osteopetrosis (marble
bone disease) results when there is incomplete bone resorption [56].

The key event that triggers the remodelling process is osteocyte programmed cell
death or apoptosis, which has been shown to depend on mechanical load [41]. Thus,
osteocyte death characterizes the regions under high loading stress that require new
bone production, as well as the zones with lower mechanical demand where bone
should be removed. In the case of micro-fractures osteocyte death is induced, which
leads to subsequent activation of the bone remodelling process. In contrast, in large
bone fractures, osteocytes release alert signals that recruit immune cells and can
lead to an inflammatory response, thus initiating an alternative pathway of bone

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 16

Figure 1.6: Coupling between bone resorption and formation. The equilibrium between both
processes is condition for the repair of microscopical skeletal damage that occurs as a result of
constant ground impact. Figure courtesy of Dr. José Manuel Lopez.

formation that involves other cell types [4].

To this day the mechanisms of cell communication within a BMU remain only
partially understood. It’s is clear however, that bone remodelling has to be tightly
regulated. For instance, an osteoblast should not start secreting osteoid matrix
before osteoclasts destruction task is finished. This raises at once the question of
how is this coordination achieved, and how could we possibly ascertain the internal
circuitry within a BMU that keeps it operative when remodelling is needed, and
shuts it off immediately afterwards. Those questions, once being clarified, would
help us to understand the coupling between bone resorption and formation. This in
turn would greatly extend our current knowledge of bone homeostasis, the strategist
to preserve it, and possible methods to restore it when necessary.

1.4 Fracture Healing

The reparation of both bone micro-fractures and larger bone fractures are complex
processes that recapitulate aspects of embriological skeletal development [26], but
in the case of larger fractures there is a particular sequence of facts that have been
well documented. It consists of three overlaping phases termed as inflammatory,
reparative and remodelling, each characterized by specific cellular, histological and
biochemical features. We next briefly sketch them.

Right after bone breaking the resulting gap is filled with blood. Thrombin gener-

Universidad Complutense de Madrid



1.4. FRACTURE HEALING 17

ated by blood vessel disruption induces the cleavage of fibrinogen present in blood to
yield fibrin monomers, which actively polymerize afterwards. In this manner a fibrin
net is generated that in 6 to 8 hours connects both sides of the fracture, and provides
a first scaffold on which bone will be eventually produced [26], [16]. Soon after bone
fracture, a number of growth factors (BMP, VEGF, PDGF, ...) are liberated within
the fracture gap [5]. Such substances will mediate subsequent steps in bone repair
(Figure 1.7, A, B).

As soon as a fibrin clot begins to be formed, new functional blood vessels as
well as Mesenchimal Stem Cells (MSC) are recruited from the fracture boundaries.
MSC move in the wake of the new vascular bed created by incoming blood vessels.
The latter can be seen within 6 to 12 hours near fracture boundaries, and between
24 to 48 hours they have invaded the whole fracture gap [50].

MSCmove inwards from the fracture boundaries, and in doing so they experience
both symmetric and asymmetric divisions [36]. Symmetric division, whereby a MSC
gives raise to two MSC, takes place whenever their available initial number is not
enough to keep the healing process on time. On the other hand, when a MSC
reaches the fibrin clot it undergoes asymmetric division, thus giving raise to one
MSC and to one differentiated cell, whose lineage depends on the density of the net
on which division takes place [10], [15], [48]. Cell replication is stimulated by some of
the growth factors present, whose action is dependent on the unfolding of the fibrin
clot (Figure 1.7 C, D). In particular, anchorage into the net significantly extends the
half-life of growth factors with respect to that of freely diffusing ones [53], [9].

Specifically, the first type of asymmetric division takes place over a low-density
fibrin clot and gives raise to fibroblasts. Once attached to the net, fibroblasts secrete
a number of substances (collagens) which stick to the clot and change its mechanical
properties [15], [48]. The resulting network (called fibrous callus) has larger me-
chanical stability than the earlier fibrin net, although not enough to ensure working
stability to the repairing bone (Figure 1.7, E, F, G).

In a period ranging from 48 to 72 hours, the fracture gap is filled with MSC and
their progenies ([32], [12]), which change as the composition of the net expanding
the fracture gap does so. In particular, upon sticking to collagen-enriched fibrin
net, MSC again undergo asymmetric division that now gives raise to one MSC and
one chondrocyte ([36], [10]). This new type of cell interacts with the existing net
and produces a more stable compound, called cartilaginous callus (Figure 1.7, H, I,
J). Bone formation now begins in a subsequent remodelling phase. More precisely,
cartilaginous callus is replaced first by trabecular and then by lamellar bone ([53])
through a process similar to endochondral ossification EO (Figure 1.7, K, L, M, N).

The micro-fracture healing is of intramembranous type and does not includes an
inflammatory phase. In contrast with large fracture, in micro-fracture repair the
new bone is formed without resorting to the formation of a fibrin-collagen scaffold
and is carried out by the action of the BMUs.

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 18

Figure 1.7: Schematic representation of the fracture healing process. It is divided into four stages:
the initial inflammatory response (A-B), soft callus formation (C-E), hard callus formation (F-I)
and bone remodeling (J-K). Figure courtesy of Dr. José Manuel Lopez

Universidad Complutense de Madrid



1.5. MATHEMATICAL MODELS 19

1.5 Mathematical models

As in any other experimental science, what we know about a phenomenon is not
the phenomenon itself, but what our understanding, mediated by our language, is
able to formulate about the problem studied [33]. In the particular cases of fracture
healing and bone remodeling , there are several conceptual models that help us to
understand such processes. Some of them will be shortly described below.

Before doing that, a few remarks are in order. To begin with, when deriving
mathematical models , two possible (not necessarily conflicting) goals may be pur-
sued. First, one might aim at proposing a set of equations whose solutions should
fit a given set of data as closely as possible. One might instead intend that solutions
to the model could explain the phenomenon, in that they could be used to make
testable predictions that could be falsified by experiments. The first approach looks
particularly appealing, and has perhaps been more actively pursued than the second
one. However, the history of science has many examples of mathematical models
that turned out to be wrong, even thoug they nicely fit data, which they did by re-
sorting to a large number of equations or to a large number of parameters appearing
in them. A particularly well known case is provided by Ptolemy’s theory of epicicles
introduced to fit planetary motions, and many other examples can be found in differ-
ent scientific fields, as for instance in theoretical ecology [27]. Such attempts, which
showed great ingenuity, often managed to match actual observations. However, they
provided no explanation about the actual reasons for the facts observed, so that no
real insight on the processes considered was gained for their use.

In this work we shall adopt a different stance, and look for model simplicity as a
basic requirement to better understand the phenomena studied. We will thus keep in
line with Ockham famous dictum : among competing hypotheses, that with the fewest
assumptions should be selected [28]. Therefore , in each of the cases addressed we
will make explicit a minimal number of assumptions (that will always be biologically
motivated) required to explain the process being dealt with. Before doing that,
however, a short review of previous models on bone remodeling and fracture repair
will be discussed below. For the ease of presentation, such models can be classified
in different types as follows.

1.5.1 Mechano-regulatory models

It is widely accepted that mathematical modeling of bone function can be traced
back to the pioneering work by Roux [47] and Wolff [55] in the XIX century. In
particular, W. Roux hypothesized that the biological processes of bone formation and
remodeling are regulated by cell signals generated by mechanical loading, a statement
considered as the fundamental premise for mechano-biology [38]. He further proposed
that cells within a tissue compete for functional stimulus. On his side, J. Wolff
further elaborated the idea that changes in bone structure and shape are a direct
consequence of changes on the mechanical stresses acting on bone, a statement that
goes (under various versions) under the generic name of Wolff’s law. Such general
postulate has not remained free of controversy. For instance, according to Cowin [11],

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 20

it considers as identical objects that are not so, as for instance stress trajectories in
a homogeneous isotropic material and the trabecullar structure of cancellous bone.
This fact notwithstanding, the assumption of a close link between stress induced and
bone structure has remained an inspiring principle in many later studies.

As a matter of fact, the literature on mathematical modeling of bone function
and repair has been steadily increasing since the sixties of last century, when the
availability of fast and efficient computers allowed for an easier comparison of nu-
merical solutions and experimental data. To classify such models, it has become
usual to consider the type of stimuli, mechanical or biological, which is considered as
more relevant in each case [25]. Such cases will be respectively termed as mechano-
regulatory and bio-regulatory.

Concerning the first type, a breaktrough is usually credited to the work done by
Pauwels about 1960 [44]. Based on Roux’s ideas, he proposed a theoretical framework
according to which mechanical forces acting on tissues determine cell differentiation
pathways on such tissues . In particular, he postulated that distortional shear stress
is a specific stimulus for the development of collagen fibers, and that hydrostatic
compressive stress is a specific stimulus for cartilage formation. According to this
author, bone formation requires of a stable, low-strain mechanical environment and
endochondral bone formation only occurs after soft tissues have stabilized their en-
vironment to that end. Pauwels theory has some limitations though. For instance,
no specific stimulus for bone formation is proposed, and his author was not able to
measure (or to estimate) in detail actual bone strains or stresses.

Further steps were provided by later models due to Carter et al. [6] and to Claes
and Heigele [7]. In these works biological tissues were considered as linear elastic
materials, and the formation of bone, cartilage and fibrous cartilage was related to
hydrostatic stress and to principal strains. In [6] the case of cyclic mechanical loading
was also considered. Different tissue differentiation mechanisms were proposed by
Prendergast et al. [45], who considered biological tissues as poroelastic materials
and linked differentiation into bone, cartilage or fibrous tissue to differences in fluid
flow properties and distortional strains. Subsequently, numerical studies based on
the approach proposed in [45] were conducted by Huiskes et al. in [30]. Finally,
a fuzzy logic theoretical model was proposed by Ament and Hofer in [1], in which
mechanical stimulus was evaluated by means of a strain energy.

A common feature of the works just mentioned is the extensive use of numerical
simulations to compare results corresponding to different bone repair situations with
experimental data.

1.5.2 Bio-regulatory models

A different approach to bone remodeling was introduced in work by Bailon-Plaza
and van der Meulen [2]. These authors proposed a mathematical model where the
effect of chemical cues, particularly growth factors, in bone remodeling was intro-
duced. Ever since, models where chemical signals are the most relevant players
retained are customarily termed as bio-regulatory. In [2], the action of growth fac-
tors on intramembranous and endochondral ossification was discussed. To do that,

Universidad Complutense de Madrid



1.5. MATHEMATICAL MODELS 21

authors formulated a PDE system involving seven variables, namely the densities
of mesenchymal stem cells (MSCs), osteoblasts and chondrocytes, concentrations
of chondrogenic and osteogenic growth factors, and densities of connective/cartilage
extracellular matrix (ECM) and bone ECM. Authors considered processes as cell dif-
ferentiation, proliferation, migration and death, as well as synthesis and resorption
(destruction) of bone tissue. A number of interesting conclusions about the regula-
tory role played by growth factors in bone healing were drawn from the numerical
simulations performed.

In 2008, Geris et al. presented in [24] a model based on previous work in [23] which
includes additional features as the impact of angiogenesis (that is, the recruitment of
new vessels from a pre-existing vasculature) on bone repair. To that end, additional
variables were introduced to those already considered in [23]: density of endothelial
cells, fibroblasts concentration, vascular matrix and fibrous tissue densities as well
as a generic angiogenic growth factor. In that way a model was arrived at containing
12 variables in a similar number of coupled PDEs with 73 parameters in all.

The cases considered in our two previous paragraphs focused mainly on mechan-
ical or biological regulation of ossification and bone repair. It is natural to consider
the possibility of integrating both approaches, thus leading to what are known as
mechano-bio-regulatory models. A significant step in this direction was provided by
Bailon-Plaza and van der Meulen in [3]. Building on their previous work [2], they
presented there a framework which specifically incorporates the effect of mechanical
stimulus on cell activity. This approach, which naturally results in increasing levels
of model complexity, has been later pursued by many authors ( see for instance[39],
[46]). Besides there are models for bone remodelling that represent remodeling on
the scale of an individual basic multicellular unit (BMU) ( see [35], [35]), and models
that extending the later ones, include the osteocyre population to acomodate the
effects of mechanical loading [40]. However, we shall refrain from adding further
details here since we intend to explore a different approach in this memoir. In-
stead, further relevant information of previous work will be provided in the following
chapters where adequate.

It follows from our previous remarks that the quest for comprehensive models,
where all biological process known to date are included , necessarily leads to large
systems of equations with an even larger number of parameters (representing reac-
tion rates, diffusion coefficients, etc) which are basically unknown. It is natural to
select them so that a particular experiment could be satisfactorily fitted. However,
most often the values thus obtained should be changed if a different experiment
is considered, and in general little insight can be gained in this manner about the
internal logic of the process under consideration. As we have already mentioned,
one could however try a different approach, and look instead for the minimal set of
biological ingredients that have to be retained in a model so that it can reproduce
what is actually observed. See for instance [19] for an example of a bio-regulatory
model in this vein. See also [29] for a survey of mathematical techniques to deal with
bio-regulatory models.

Universidad Complutense de Madrid



CHAPTER 1. PRELIMINARIES 22

Universidad Complutense de Madrid



Bibliography

[1] Ament, C. & Hofer, E. P. 2000 A fuzzy logic model of fracture healing, Journal
of Biomechanics, 33, (2000), 961-968. (doi:10.1016/S0021-9290(00)00049-X)

[2] Bailon-Plaza, A. & van der Meulen, M. C. H. 2001 A Mathematical framework
to study the efects of growth factor influences on fracture healing, Journal of
Theoretical Biology., 212, (2001), 191-209. (doi:10.1006/jtbi.2001.2372)

[3] Bailon-Plaza, A. & van der Meulen, M. C. H. 2003 Beneficial effects of mod-
erate, early loading and adverse effects of delayed or excessive loading on bone
healing, Journal of Biomechanics., 36, (2003), 1069-1077. (doi:10.1016/S0021-
9290(03)00117-9)

[4] J.P. Bidwell, J. Yang & A. G. Robling 2008. Is HMGB1 an osteocyte alarmin?.
J. Cell. Biochem. 103/6 1671-1680.

[5] Bolander, M.E. 1992. Regulation of fracture repair by growth factors. Proc Soc
Exp Biol Med. 200, 165-170.

[6] Carter, D. R., Blenman, P. R. & Beaupre, G. S. 1988 Correlations between me-
chanical stress history and tissue differentiation in initial fracture healing., J.
Orthop. Res., 6, 736-748. (doi:10.1002/jor.1100060517)

[7] Claes, L., Wilke, H.-J., Augat, P., Rubenacker, S., & Margevicius, K., 1995 Efect
of dynamization of gap healing of diaphyseal fractures under external fixation,
Clinical Biomechanics ., 8, (1995), 227-234.

[8] Claes, L. E. & Heigele, C. A. 1999Magnitudes of local stress and strain along bony
surfaces predict the course and type of fracture healing, Journal of Biomechanics.,
32, (1999), 255-266. (doi:10.1016/S00219290(98)00153-5)

[9] Cole, A.A., & Walters, L.M. 1987. Tartrate-resistant acid phosphatase in bone
and cartilage following decalcification and cold-embedding in plastic. J. His-
tochem. Cytochem. 35, 203-206.

23



BIBLIOGRAPHY 24

[10] Colnot, C., Huang, S., Helms, J. 2006. Analyzing the cellular contribution of
bone marrow to fracture healing using bone marrow transplantation in mice.
Biochem. Biophys Res Commun, 350, 557-561

[11] Cowin, S. 2001. The false premise of Wollf ’s law. In: Cowin, S.C. (Ed.), Bone
Biomechanics Handbook. CRC-Press, Boca Raton, pp. 1?14 (Chapter 30).

[12] Deckers, M.M.L., van Bezooijen, R.L., van der Horst, G., Hoogendam, J., van
der Bent, C., Papapoulos, S.E. et al. 2002. Bone morphogenetic proteins stimu-
late angiogenesis through osteoblast-derived vascular endothelial growth factor A.
Endocrinology, 143, 1545-1553.

[13] Doll, B. 2005 Develpmental Biology in the Skeleton system. In Bone Tissue
Engineering. Edited by Hollinger, J., Einhorn, T. Doll, B. & Sfeir, C. CRC Press.

[14] Donahue, H. J., Siedlecky, C. A. & Vogler, E. 2005 Osteoblastic and Osteo-
cytic Biology in Bone Tissue Engineering. In Bone Tissue Engineering. Edited
by Hollinger, J., Einhorn, T. Doll, B. & Sfeir, C. CRC Press.

[15] Eghbali-Fatourechi, J., Lamsam, D., Fraser, D., Nagel, B.L., Riggs, S. & Khosla
2005. Circulating osteoblast-lineage cells in humans. N Engl J Med. 352, 1959-
1966.

[16] Einhorn, T.A. 1998. The cell and molecular biology of fracture Healing.
Clin.Orthop.Relat. Res 355, S7-21.

[17] Encyclopaedia Britannica (Online). Documenting electronic sources
on the internet. (Date of consult: october 2014). Available in:
http://global.britannica.com/EBchecked/topic/434325/osteon

[18] Erik Fink Eriksen 2010 Cellular mechanisms of bone remodelling Rev Endocr
Metab Disor 11 219-227

[19] Fasano, A. Herrero, M., Lopez, J.M., & Medina, E. 2010 On the dynamics of
the growth plate in primary ossification. J. Theor. Biol. 265 543-553.

[20] Ferguson, C., Miclau, T., Alpern, E. & Helms, J. 1999 Does adult fracture repair
recapulate embrionic skeletal formation? Mech. Dev. 87 57-66

[21] Frost H. M. 1964 Dynamics of bone remodelling. Bone Biodyn. 315-333. (Book
Chapter)

[22] Gay, C.V. 2005 The osteoclast In Bone Tissue Engineering. Edited by Hollinger,
J., Einhorn, T. Doll, B. & Sfeir, C. CRC Press.

[23] Geris, L. Gerisch, A., Maes, C., Carmeliet, G., Weiner, R., Van Oosterwyck,
H., & Vander Sloten. 2006 Mathematical model of fracture healing in mice: com-
parision between experimental data and numerical simulations. Med Biol Eng
Comput 44, 280-289. (doi:10.1007/s11517-006-0040-6)

Universidad Complutense de Madrid



BIBLIOGRAPHY 25

[24] Geris, L., Gerisch, A., Vander Sloten, J., Weiner, R. & Van Oosterwyck, H.
2008 Angiogenesis in bone fracture healing: a bioregulatory model. Journal of
Theoretical Biology 251, 137-158. (doi:10.1016/j.jtbi.2007.11.008)

[25] Geris, L., Vander Sloten, J., & Van Oosterwyck, H. 2009 In silico biology of bone
modelling and remodelling, Philosophical Transactions of the Royal Society., 367,
(2009), 2031-2053. (doi:10.1098/rsta.2008.0293)

[26] Gerstenfeld, L.C., Cullinane, D.M., Barnes, G.L., Graves, D.T., & Einhorn,
T.A. 2003.Fracture healing as a post-natal developmental process: molecular, spa-
tial, and temporal aspects of its regulation. J Cell Biochem. 88, 873-884.

[27] Ginzburg, L. R. & Jensen, C. X. J. Rules of thumb for judging ecological theories.
Trends in Ecology and Evolution. 19/3 121-126

[28] Hamilton, Sir William. 1860. Lectures on Metaphysics and Logic. Oxford Uni-
versity Press.

[29] Herrero, M. & Lopez, J.M. 2005 Bone formation: biological aspects and modeling
problems. Theor. Med. 6/1 41-55.

[30] Huiskes,R., van Driel,W.D., Prendergast,P.J.,& Soballe,K. 1997. A biomechani-
cal regulatory model for periprosthetic fibrous-tissue differentiation. J. Mater. Sci.
Mater. Med. 8 785-788. (doi:10.1023/A:1018520914512)

[31] Huiskes, R. 2000 If bone is the answer, what is the question?. J. Anat. 197.
145-156.

[32] Kalfas, I.H. 2001. Principles of bone healing. Neurosurg Focus, 10.

[33] Kant, I. 1788. Critique of Pure Reason. Translated into english by Norman
Kemp Smith 1929. Cambridge University Press.

[34] Komarova, S., Smith, R., Dixon, S., Sims, S., & Wall, L. 2003 Mathematical
model predicts a critical role for osteoclast autocrine regulation in the control of
bone remodelling. Bone 33 225-234.

[35] Komarova, S. 2005 Mathematical model of paracrine interactions between osteo-
clasts and osteoblasts predicts anabolic action on parathyroid hormone on bone.
J. Endocrinol. 146 3589-3595

[36] Malizos, N. & Papatheodorou, L.K. 2005. The healing potential of the perios-
teum. Molecular aspects. Injury 36 (Suppl 3), S13-S19.

[37] Marks, S. C. & Odgren, P. R. 2002. Structure and devopment of the skeleton.
In Principles of Bone Biology 2 ed. Edited by Bilezikian, J. Raisz, L. & Rodan,
G. Academic Press.

[38] van der Meulen, M.C., & Huiskes, R., 2002 Why mechanobiology ? A survey
article J. Biomech. 35 401-414.

Universidad Complutense de Madrid



BIBLIOGRAPHY 26

[39] Moreo, P., Garcia Aznar, J., Doblare, M. 2009. Bone ingrowth on the surface
of endosseous implants. part i: mathematical model. J. Theor. Biol 260, 1–12
(2009)

[40] Moroz, A., Crane, M., Smith, G., & Wimpenny, D. Phenomenological model
of bone remodelling cycle containing osteocytes regulation loop. Biosystems. 84
183-190.

[41] B.S. Noble, N. Peet, H.Y. Stevens, A. Brabbs et al. 2003. Mechanical loading
: biphasic osteocyte survival and targeting of osteoclasts for bone destruction in
rat cortical bone. Am. J. Phys-Cell Phys. 284/4 C934-C943.

[42] Overgaard, S. 2000 Calcium phosphatae coatings for fixation of bone im-
plants., Acta Orthopaedica Escandinavica Suplementum, 297, (2000), 1-74.
(doi:10.1080/000164700753759574)

[43] Parfitt, A. M. 2002 Targeted and not targeted bone remodelling: Relationship to
basic multicellular unit origination and progression. Bone 30/1, 5-7.

[44] Pauwels, F. 1960 Eine neue Theorie uber den Einflus mechanischer Reize auf die
Differenzierung der Stutzgewebe. Z. Anat. Entwicklungsgeschichte 121, 478-515.
(doi:10.1007/BF00523401)

[45] Prendergast, P. J., Huiskes, R. & Soballe, K. 1997 Biophysical stimuli on
cells during tissue differentiation at implant interfaces., J. Biomech. 30 539-548.
(doi:10.1016/S0021-9290(96)00140-6)

[46] Prokharau, P. Vermolen, F. & Garcia-Aznar, J. M. 2012 Model for direct bone
apposition on pre-existing surfaces, during peri-implant osseointegration., Journal
of Theoretical Biology., 304, (2012), 131-142.

[47] Roux, W. 1881 Der Kampf der Teile im Organismus: ein Beitrag zur Vervoll-
standingung der mechanischen Zweckmassigkeitslehre. Engelmann, Leipzig.

[48] Rumi, M.N., Deol, G.S., Singapuri, K.P., & Pellegrini Jr. V.D. 2005. The ori-
gin of osteo- progenitor cells responsible for heterotopic ossification following hip
surgery: an animal model in the Rabbit. J Orthop Res. 23, 34-40.

[49] Sims, N. & Martin, T. J. 2014 Coupling the activities of bone formation and
resorption: a multitude of signals within the basic multicellular unit. BoneKEy
Reports 3, Article Number 481. 1-10.

[50] Shapiro, F. 2008 Bone development and its relation to fracture repair: The role
of mesenchymal osteoblasts and surface osteoblasts. European Cells and Materi-
als. 15 56-73.

[51] Steiniche, T., & Hauge, E. M. 2003 Normal structure and function of bone. In
Handbook of Histology for Bone and cartilage. Humana Press. New Jersey.

Universidad Complutense de Madrid



BIBLIOGRAPHY 27

[52] Strom, O. Borgstrom, F., Kanis, J. A. Compston, J. Cooper, C. McCloskey, E.
V. & Johnsson, B. 2011. Osteoporosis: burden, health care provision and opportu-
nities in the European Union. Arch Osteoporos. 6. 59-155 DOI 10.1007/s11657-
011-0060-1.

[53] Takahara, M., Naruse, T., Takagi, M., Orui, H., & Ogino, T. 2004. Matrix
metalloproteinase-9 expression, tartrate-resistant acid phosphatase activity, and
DNA fragmentation in vascular and cellular invasion into cartilage preceding
primary endochondral ossification in long bones. J Orthop Res. 22, 1050-1057.

[54] U.S. Department of Health and Human Services 2004. Bone Health
and osteoporosis: A Report of the Surgeon General. Available at:
http://www.ncbi.nlm.nih.gov/books/NBK45513/ .

[55] Wolf, J. 1892. Das Gesetz der Transformation der Knochen. Verlag von August.
Berlin. Translated as: The Law of bone remodelling. Springer Verlag. Berlin, 1986.

[56] Zaidi, M. 2007 Skeletal remodelling in health and disease. Nature Medicine. 13/7
791-801

Universidad Complutense de Madrid



BIBLIOGRAPHY 28

Universidad Complutense de Madrid



Chapter 2
Bone remodelling

2.1 Introduction

It is well known that bone microscopic structure continuously adapt to meet the
resulting mechanical requirements. This adaptation is achieved by means of a bone
remodeling process that is maintained throughout life, and whose control relies on
osteocytes, the most abundant cell type in bone [1, 2].

Osteocytes act as mechanosensors that monitor mechanical stress within bone
tissues, reacting to changes in both the amount and the direction of loading on
bone and ensuring that remodeling only takes place where needed. The key event
that triggers the remodeling process is osteocyte programmed cell death or apop-
tosis, which has been shown to depend on mechanical load. The relation between
osteocyte apoptosis and load applied is U-shaped, i.e., mechanical stresses within a
normal physiological range prevent it, whereas those above or below this range favor
it [3, 4, 5]. Thus, osteocyte apoptosis is a marker both for regions under high loading
stress that require new bone production, and for zones with lower mechanical de-
mand where bone should be removed. Micro-fractures, which are continuously being
produced as a consequence of physical exercise in healthy individuals, also result in
osteocyte apoptosis and subsequent activation of the bone remodeling process. In
bone fractures, a considerable amount of osteocytes are compromised and alert sig-
nals are produced that recruit immune cells and lead to an inflammatory response.
In this manner an alternative pathway of bone formation is started where other cell
types are involved [6]. We shall not deal with this second situation here, since we
will be concerned with homeostatic bone remodeling at comparatively smaller scales.
The way in which this process is carried out is now discussed below.

Once triggered by osteocyte apoptosis, a bone remodeling process is executed
by organic teams called Basic Multicellular Units (BMU) [7]. Each BMU consists
of several morphologically and functionally different cell types (mainly osteocytes,
osteoblasts and osteoclasts) that act co-ordinately in small cavities, called Bone re-

29



CHAPTER 2. BONE REMODELLING 30

modeling Compartments (BRC) where remodeling takes place [8, 9]. The behaviours
of the different cell types involved can be summarized as follows. A normal pres-
ence of osteocytes keeps osteoblasts deactivated and precludes osteoclast precursors
from being recruited. When a group of osteocytes undergoes apoptosis, a drop in
inhibitory action results that has several consequences. In particular, it leads to
osteoblasts activation and to the attraction of osteoclast precursors arising from the
bone marrow. These precursors subsequently differentiate into mature osteoclasts,
which start destroying surrounding bone, leading to the appearance of the so-called
cutting cones [10]. In a later stage of the process, activated osteoblasts, marching
in the wake of bone-destroying osteoclasts, will replenish the cavity left behind the
latter with an osteoid matrix. Some of the active osteoblasts remain trapped in the
bone matrix and differentiate into osteocytes [11], which in turn will be instrumental
in the subsequent mineralization of the matrix surrounding them, thus concluding
local bone remodeling. The whole process is summarized in Figure 1.

Figure 2.1: A sketch of BMU operation after a group of osteocytes has
undergone apoptosis near bone surface. Left) Bone interior (in gray) contains
a number of functional osteocytes (OCY, in green). Apoptotic osteocytes (aOCY,
in black) cease to release osteoblasts (OBL) inhibitors, which in particular leads to
the activation of OBL placed at the bone surface. Center) A remodeling process
starts. As a further consequence of the drop in OCY-induced inhibitory effects,
Osteoclast precursors are recruited and differentiate into active osteoclasts (OCL,
in red). OCL then form cutting cones that start the destruction (resorption) of
surrounding bone. Later on, activated osteoblasts (OBA) move inwards in the trail
open by OCL. Right). When the bone-remodeling compartment (BRC) reaches a
large enough size, OCL die by apoptosis and OBA secrete osteoid matrix. OBA
then differentiate into OCY and remain trapped in that matrix, which is eventually
mineralized, thus completing the formation of new bone.

The process just recalled is exquisitely regulated. For instance, in addition to
taking place only where is needed, the amount of new bone generated should equal
that of old bone destroyed. Achieving such a balance requires of a coordinated action
of all cell types involved. To this day, and in spite of substantial progress achieved

Universidad Complutense de Madrid



2.1. INTRODUCTION 31

during the last 50 years, the mechanisms of cell communication within a BMU remain
only partially understood. It is known however that such regulation is of a local
nature since remodeling, which is continuously taking place, is simultaneously and
independently occurring in different parts of the skeleton [12]. The process under
consideration is a consequence of individual cell decisions. In fact, at any place where
the process is going on, individual cells within a BMU are not aware of the pace at
which remodeling is proceeding. Nevertheless, to efficiently fulfill its function, any
such cell should perform its task according to a precise schedule. For instance, an
osteoblast should not start secreting osteoid matrix before osteoclasts destruction
task is finished. This naturally raises the question of how is this timing achieved,
and what internal circuitry within a BMU keeps such unit operative when remodeling
is needed, and shuts it off immediately afterwards. We shall address these issues in
the following sections.

Specifically, in this work we formulate and analyze a simple, space-dependent
mathematical model that is able to describe BMU operation. Simplicity is here a
key issue. We assume that any cell type within a BMU can only take a very limited
number of actions (divide, die, differentiate) each being blocked within that cell by a
specific inhibitor molecule. The internal dynamics of such inhibitors is coupled with
the external medium via a minimal number of cell-to-cell signaling cues, which upon
binding with membrane receptors act as sources or sinks for the internal inhibitors.
The first of them which falls below a threshold value determines a decision within the
previous list in each cell. We will show that only three signaling molecules, acting on
cells that can select one between at most three options, are enough to keep a BMU
fully operational. Other than choosing a given cell fate over another, cells in a BMU
move in response to environmental signals. For instance, osteoclasts form cutting
cones by moving inwards into the bone to be resorpted. On the other hand, active
osteoblasts lying behind cutting cones will eventually move in the opposite direction
during bone formation stage. This dynamic aspect of BMU operation will also be
included in our mathematical model.

Once our approach has been stated, a few remarks are in order. To begin with,
one may wonder how three signaling cues can possibly be enough for our purposes,
taking into account that a much larger set of mediating molecules have been de-
scribed in the context of bone remodeling (some of them will be recalled below for
convenience of the reader). We believe that this is related to some features of the
signaling network involved that are far from being fully understood as yet. Namely,
we refer to multi-functionality (the fact that a given signal may induce quite different
effects in different cell types ) and redundancy ( the presence of several molecules
that may produce the same effect in a given cell type). We do not aim here at pro-
viding a priority list between the signals described. What we show instead is that
three such cues are enough to regulate the internal inhibitors of the cells considered,
so that BMU operation results as an emergent property of the individual decisions
taken by any cell involved. Interestingly, such inhibitors have been described in the
literature, and the actions induced by the signals proposed have been associated to
various of the molecules identified in bone remodeling. We shall provide additional

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 32

details on these topics in our next section.
A consequence of the previous remarks is that we do not aim at providing a

comprehensive model of bone remodeling where most (if not all) possibly involved
processes are retained, but rather at identifying the minimal BMU software required
to keep such unit operative. We believe that this may be a useful complement to
previous mathematical models available (see for instance [13, 14, 15] and references
therein). In these works, great attention is paid to include as many biological data
as possible in the models considered, which in particular leads to difficult problems
related to actual parameter fitting for the processes involved. In spite of such obsta-
cles, and as a consequence of authors ingenuity, quite interesting results have been
obtained in such (and other) references. We hope that a combined use of both ap-
proaches (comprehensive and minimality -driven) could shed further light into our
understanding of bone remodeling , and particularly on its control when such type
of action would be required .

We conclude this introduction by describing the plan of this work. We first
recall in Section 2.2 below some known results about cell types involved in BMU
functioning, their internal inhibitors and some of the main chemical signals that
have been described to mediate bone remodeling. In Section 2.3 we then formulate
our mathematical model which is subsequently simulated in Section 2.4 to yield a
number of results. In particular, we discuss there the manner in which the number
of cells in a BMU is determined by the BRC size, describe how osteoclast lifespan
depends on the depth of the BRC and show that the BMU circuit proposed is robust
to external noise. The paper then concludes with a Section 2.5, where our results
are discussed and possible future directions are shortly described.

2.2 A minimal software for BMU operation

In this section we shall briefly recall some of the chemical signals that are known
to be involved in bone remodeling by BMUs. Specifically, a list will be provided of
the main molecular mediators in BMU operation currently known, and of the effects
they induce in the different cell types involved. In-depth reviews can be found, for
instance, in [8, 16, 17, 12].

2.2.1 The role of molecular signals in bone remodeling

A considerable number of chemical signals involved in BMU actions have already
been identified. For instance, sclerostin is secreted by normal osteocytes and is known
to inhibit both osteoblast activation and osteoclast recruitment from precursor cells
in the bone marrow [8, 18]. This inhibitory effect prevents the activation of BMUs in
regions where bone remodeling is not needed. In fact, only when a sufficient number
of osteocytes undergo apoptosis in a given area, due for instance to micro-fractures
or changes in mechanical load, will the inhibitory effect of sclerostin be removed, thus
entailing the activation of a BMU and the subsequent initiation of a bone remodeling
event.

Universidad Complutense de Madrid



2.2. A MINIMAL SOFTWARE FOR BMU OPERATION 33

A set of signals known under the generic name of Bone Morphogenetic Proteins
(BMPs) have been described to play a key role during the first stages of bone re-
modeling. Some members of this family, such as TGF-β and TGF-α, are known to
foster differentiation of mesenchymal stem cells to osteoblasts [19, 20]. TGF-β can
also induce migration of osteoblasts to the sites of bone formation during remodeling
[21] and inhibits osteoblast apoptosis [4]. TGF-β is mainly produced by osteocytes
[18, 22], but it is also present in the bone matrix [19] and in platelets [23]. This
last factor, together with other signals such as HMGB-1 [6] provide a link with in-
flammatory processes occurring at early stages of larger fractures repair. Cytokines
such as IGF-1, released from bone matrix, seem to play a similar role in activating
osteoblast differentiation [21] and are necessary for their survival in vitro [24].

The next stage in the process consists in the recruitment of osteoclasts to the site
of bone remodeling. The main signals involved in this step are M-CSF and RANKL
[25, 26, 17], that promote the differentiation of osteoclast precursors and the survival
of activated osteoclasts. Since RANKL is mostly produced by activated osteoblasts
[8, 25], osteoclasts will only be recruited to sites were bone remodeling has already
been triggered. Furthermore, in case of a contingent activation of osteoclast precur-
sors where bone remodeling is not necessary, osteoclasts will die by apoptosis in the
absence of activated osteoblasts nearby. This provides a security mechanism that
avoids the presence of active osteoclasts out of the context of BMU operation. The
release by inactive osteoblasts and osteocytes of osteoprotegerin (OPG), a molecule
that inhibits the action of RANKL on osteoclasts [25, 26, 27], can also be understood
as a checkpoint mechanism that ensures that osteoclasts will only be activated and
kept alive in the presence of a sufficient number of activated osteoblasts. Besides
inhibiting osteoblast activation, sclerostin also induces apoptosis of active osteoclasts
[8, 18]. It may thus act as a signal that stops bone resorption when the cutting edge
has reached the required depth in each remodeling event. Finally, osteoid matrix pro-
duction by active osteoblasts, as well as differentiation of osteoblasts into osteocytes
seem to be cell-density dependent [28], and have been suggested to be mediated by
connexin, a molecule that circulates through gap junction communications between
osteoblasts, as well as by sclerostin [16, 29].

Concerning signalling effects, one also has to bear in mind that: I) different signals
can have redundant effects; for example both TGF-β and IGF-1 induce osteoblast
differentiation [21], II) a same signal can have different, even opposite, effects on
different cell types. For instance TGF-β, FGF and PDGF activate osteoblast and
inhibit osteoclast action [23] and G-CSF is known to induce apoptosis and to inhibit
differentiation in osteoblasts [26], and III) signals are not always provided by chem-
icals. In this context we have already remarked that osteocytes act as sensors that
respond to changes in mechanical stress in bone [27, 31].

The list just provided of molecular mediators and their effects on BMU cell
types, which is summarized in Table 1 below, is far from being complete. Moreover,
knowing the identity of molecular mediators or describing their effects in qualitative
terms is not enough to explain BMU operation during bone remodeling. In order
to understand how a coherent collective plan of action emerges at a multicellular

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 34

scale, it is also necessary to take into account quantitative aspects of the process.
Indeed, for any given set of signals involved, the amount of bone to be resorpted
and produced in different remodeling events can be highly variable [8]. This implies
in turn that the number and activity of cells recruited in a BMU should change to
suit the needs of each particular remodeling process [22, 23, 30, 32]. Notice also that
the references quoted do not focus on quantitative aspects of the process, such as
the rate of the production of signals by BMU cell types, or their effect at different
concentrations.

Signal Source OCL OBL OBA

HMGB-1 aOCY
Recruitment
Chemotaxis

Recruitment
Chemotaxis

Chemotaxis

RANK, OPG
M-CSF

OBA
aOCY

Recruitment
Activation
Survival

- -

Sclerostin OCY Inhibition
No proliferation
No differentiation

Apoptosis

No proliferation
No differentiation

Apoptosis
TGF-β
BMPs
PDGF
IGF
FGF

OCY
OBL
OBA

Inhibition

Chemotaxis
Recruitment
Differentiation
Proliferation
Survival

Chemotaxis
Differentiation
Proliferation
Survival

Table 2.1: A brief summary of the main signals involved in BMU operation as
described in the literature (see references in the text)

In our view, the key to understand BMU operation is signal integration by indi-
vidual cells, as discussed in the next section.

2.2.2 Signal integration by individual cells. Inhibitory pro-
teins

We now postulate our basic modelling assumptions, which can be summarized as
follows. Cells within a BMU integrate the molecular signals in their immediate
surroundings and, based on the information conveyed by these signals, they chose one
among a very limited set of actions, namely differentiation, cell division, migration
and apoptosis. In addition, we postulate that inhibitory proteins, that block the
progression of these actions within each cell type, mediate cell decisions in any bone-
remodeling event.

This last assumption is illustrated by the well-known behaviour of two inhibitory
proteins, Rb and Bcl-2. The retinoblastoma protein (Rb) arrests progression into the
cell cycle, whereas the B-cell lymphome-2 protein (Bcl-2) precludes the initiation of
the apoptosis program in most cell types, including osteocytes and other stromal cells
[27]. More precisely, Rb binds to transcription factors of the E2F family, preventing

Universidad Complutense de Madrid



2.3. MATHEMATICAL MODELLING 35

the progression of the cell cycle to the synthesis stage [33]. When the amount of active
Rb falls below a critical threshold, a no-return point is reached (the Restriction
Point of the cell cycle) that irreversibly leads to cell division. Analogously, Bcl-2
precludes the release of cytochrome c through the mitochondrial outer membrane,
thus avoiding the initiation of the apoptosis program. If Bcl-2 is depleted beneath a
suitable value, its inhibitory action is lost and the cell is committed to die [34, 23].
Hence a competition between inhibitory molecules results in a mechanism of cell fate
control: the first of these inhibitors (Rb or Bcl-2) that falls below its corresponding
threshold value determines the cell decision to divide or die and, importantly, the
timing of such decision (see Figure 2.A). This mechanism also provides an explicit
link between external signals and cell decisions, given that membrane receptors are
known to modulate the evolution of inhibitory molecules within the cell [34]. As a
matter of fact, it has been recently shown that the interplay between receptor/signal
interaction and the internal dynamics of Rb and Bcl-2 suffices to explain the onset
of emergent population properties (as clonal expansion and contraction) in the case
of immune response to acute infections [35].

We suggest that a similar approach can be implemented to unravel a BMU op-
erational mechanism. For instance, the existence of inhibitory proteins controlling
cellular processes during bone remodeling is well documented in the literature. In
fact, two restriction points have been identified in the differentiation program of os-
teoblasts (marking respectively the transition to activated osteoblast and osteocyte
types) that lead to the arrest of the process in the absence of the right signals [23].
One of these restriction points is determined by the inhibition of the transcription
factor Runx2 [36] by the action of proteins of the Twist family [37]. Moreover,
TGF-α and BMP signalling are known to mediate Runx2 dynamics, thus modulat-
ing osteoblast differentiation [38, 39, 40]. The increased complexity derived from the
presence of more than one cell type, together with the existence of several cell fates
(division, apoptosis or alternative differentiation programs) now introduces new pos-
sibilities with respect to the basic dichotomy cell division vs. cell death considered
in [35]. However, the underlying logic can be extended to account for these new
alternatives. For instance, several inhibitory molecules can simultaneously block the
progression of alternative differentiation programs. In this case, the first inhibitor to
reach its critical threshold will unambiguously determine the fate of the cell (see Fig-
ure 2.B).We postulate that cell choices thus determined are mutually exclusive. This
seems to be the case for BMU cells. For instance, osteoblasts that start the differ-
entiation program or decide to secrete osteoid matrix do not complete the cell cycle,
and therefore do not proliferate [23]. We also remark that this mechanism allows for
one signal to induce alternative cell decisions depending on its concentration.

2.3 Mathematical modelling

Bearing in mind the multiplicity of signals recalled in Table 2.1, as well as the redun-
dancy often observed in their functioning, we propose that the effective operation
of a BMU can be modelled by means of a reduced version of the complex signalling

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 36

Figure 2.2: A mechanism of cell fate determination. A) Inhibitory proteins Rb
and Bcl-2 provide a fate decision mechanism in several cell types. The first of these
molecules to reach its critical threshold determines if the cell will divide or undergo
apoptosis. For convenience all inhibitory thresholds are set equal to zero. B) The
presence of inhibitors blocking alternative differentiation programs allows to increase
the complexity of this cell fate-decision mechanism. B Left) If differentiation is
assumed to be blocked by Inhibitor 1, which happens to vanish before Rb and Bcl-
2 do, then the cell will not divide or die, but will instead undergo differentiation
program 1. B Right) Two alternative differentiation programs can be controlled by
two different inhibitors (1 and 2). If one of them (inhibitor 2 in this case) disappears
faster than the remaining molecules, the corresponding differentiation program is
selected.

network previously described. Specifically, we propose that three cell-released sig-
nals, denoted by R, S and T suffice for that purpose. The effect induced by such
signals in the cell types involved is described in Table 2.2 We will next describe the
main ingredients of the model proposed in this work.

2.3.1 Cell algorithms

In agreement with our main assumption in Section 2.2, we assume that progression
within the cell cycle, apoptosis and the initiation of differentiation programs are
initially blocked by specific inhibitory molecules. This would allow for newly formed
cells to remain in the first stage of the cell cycle (G1) before choosing a given cell fate.
During this stage, membrane receptors interacting with external signals dictate the
dynamics of inhibitors, thus controlling the eventual decision of the cell (see Figure
3). The effect of the complexes formed by membrane receptors and external cues
in the inhibition (resp. activation) of any cell fate choice takes place through an
increase (resp. decrease) in the amount of the inhibitory molecule controlling this

Universidad Complutense de Madrid



2.3. MATHEMATICAL MODELLING 37

Signal Source OCL OBL OBA

R OBA
Activation
Survival

- -

S OCY
Apoptosis
Inhibition

No proliferation
No differentiation

Apoptosis

No proliferation
No differentiation

Apoptosis

T
OCY
OBL
OBA

-
Proliferation
Differentiation

Survival
No proliferation

Table 2.2: Generic signals to be considered in our model and a list of the
actions they induce on cell types of a BMU

choice. Cell fate is decided when the first of the inhibitors considered reaches its
critical threshold.

Figure 2.3: Cell fate decision algorithm. A) The first stage of the cell cycle
is a decision phase (dec). During this phase, membrane receptors interact with
specific external signals (S1, S2 and S3 in the picture). As a result, changes in
the concentrations of inhibitory molecules blocking progression into the cell cycle
(div), apoptosis (apo) and differentiation (diff) are induced. In the particular case
considered, signal S1 competes with S2 to induce and inhibit cell division respectively.
Signal S2 induces the initiation of the apoptosis program, while signal S3 leads to
cell differentiation. Each inhibitor has a critical threshold (dashed line) that has to
be reached to trigger the corresponding cell action. B) Diagram of the cell decision
algorithm described in A). We remark that signals controlling cell fate do not need
to be produced by the cell itself. In this case, for instance, the cell releases signals
S1 and S3, while signal S2 is produced by other cell types.

The inhibitor dynamics will be assumed to be of the following form:

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 38


c′(t) = fc(R,S, T )

a′(t) = fa(R,S, T )

d′(t) = fd(R,S, T ),

(2.1)

where c(t), a(t) and d(t) respectively denote the concentrations of division, apop-
tosis and differentiation inhibitors at time t and fc, fa, and fd are functions of three
external signals S, R and T . For simplicity, we shall assume in the sequel such func-
tions to be piecewise linear, and will next describe the main details of the decision
algorithm for any cell type in a BMU.

Osteoblasts (OBL)

Osteoblasts (OBL) are initially located at the interface between bone matrix and
bone marrow. In normal conditions, osteoblasts homeostasis is maintained by a
continuous cell turnover, involving both division and apoptosis [41]. If a remodeling
process starts, osteoblasts can also differentiate into active osteoblast (OBA). OBLs
can therefore choose among three alternative programs. In this work we will not
consider OBL division and apoptosis, which are therefore assumed to balance each
other during the process, and will focus on OBL activation during bone remodeling.
Activation is known to be mainly blocked by the release of sclerostin by osteocytes
[5, 8]. This will be modeled by assuming that signal S, produced by osteoctyes,
increases the amount of a differentiation inhibitor in OBLs, denoted by dB , according
to the following dynamics:

{
d′B(x, y; t) = αS1

B S(x, y; t)− αS2

B , if dB < DB and box at (x,y) is adjacent to bone

d′B(x, y; t) = 0 otherwise,

(2.2)
where (x, y) denotes the position of the OBL, DB is the maximum amount of in-
hibitor that can accumulate within a OBL, S(x, y; t) is the amount of signal S at
location (x, y) at time t and α1

B and α2
B are positive parameters. During a bone re-

modeling process, osteoblasts remain in their original positions, delimiting the BRC
[9]. Accordingly, OBLs do not move in our model.

Active osteoblasts (OBA)

Once differentiated, OBLs become active osteoblasts (OBA). OBAs have three pos-
sible choices: cell division, cell death and differentiation into osteocytes (OCY). We
postulate that OBA division, apoptosis and differentiation are inhibited by molecules
whose concentration are denoted by cA, aA and dA. Furthermore, we assume that
following activation OBAs can only divide, which leads to the layer of OBAs to pro-
gressively cover the surface of the cutting cone (see Section 2.3.2 below). We assume
that during this process the existence of OBA that are not in contact with the bone
surface blocks both the differentiation of OBAs into OCY and their apoptosis. This

Universidad Complutense de Madrid



2.3. MATHEMATICAL MODELLING 39

can be achieved, for instance, by the interchange of molecules through gap junction
communications between cells. When all OBAs are attached to the bone matrix,
OBAs start to synthesize osteoid matrix, and can also differentiate into OCY or die,
depending on the concentration of signal S in the environment. Differentiation into
OCY will also be assumed to be blocked by cell contact with existing OCY, so that
only OBAs whose distance to OCY is larger that a threshold can differentiate. These
behaviors are modeled by means of the following equations:

{
c′A(x, y; t) = −αT

AT (x, y; t), if box at (x,y) is not adjacent to bone

c′A(x, y; t) = 0 otherwise,
(2.3)

where (x, y) denotes the position of the OBA, T (x, y; t) is the amount of signal T at
location (x, y) at time t and α1

A is a positive parameters.
When all OBAs are attached to the bone surface, cell division is no longer possible

for lack of space, and osteoid production, apoptosis and differentiation into OCY can
take place, according to these equations:


a′A(x, y; t) = −αS1

A S(x, y; t)

d′A(x, y; t) =

{
αS2

A S(x, y; t)− αS3

A , if distance to boxes occupied by OCY is less than δA

0 otherwise

(2.4)
where S(x, y; t) is the amount of signal S at location (x, y) at time t and αi

A (for i =
S1, S2, S3) are structural positive parameters and δA is a positive variable parameter.
Parameter δA allows to model the effect of cell-to-cell contact inhibition by OCYs
on OBAs activation, resulting from gap junction communications between cells [8].

Osteoclast precursors (OCP)

OBA can recruit OCP to the bone remodeling zone by secreting signal R, while signal
S inhibits OCP recruitment. We will not consider in the model OCP division and
apoptosis, and assume instead that enough OCP are available in the bone marrow
(in particular, in the locations adjacent to those occupied by OBAs). OCPs attach
to the bone surface upon their recruitment by OBAs, after which they become active
osteoclasts (OCL). We assume that OCP activation is blocked by inhibitor dP , whose
dynamics is modeled as follows:

d′P (x, y; t) = −αR
PR(x, y; t) + αS

PS(x, y; t) (2.5)

where R(x, y; t) and S(x, y; t) are respectively the amounts of signal R and S at
location (x, y) at time t and α1

P and α2
P are positive parameters.

Osteoclasts (OCL)

Once recruited to the BRC, OCPs become active osteoclasts (OCL) and start dis-
solving the bone matrix, thus leading to the formation of a cutting cone. OCLs are

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 40

assumed to be endowed with two alternative apoptosis pathways. One is controlled
by signal R and ensures that OCLs die in the absence of activated OBA. A second
pathway is controlled by signal S, produced by OCYs, and determines the depth
that will be reached by the cutting cone (the progression of OCLs into old bone is
described in Section 2.3.1). In particular, OCLs in locations with high concentrations
of S will stop digging into the bone matrix and die. We will model the apoptosis
pathway controlled by signal S by assuming that cell death is blocked by inhibitor
aC , whose dynamics is given by following equations:{

a′C(x, y; t) = αS1

C if S(x, y; t) > αS2

C

a′C(x, y; t) = 0 otherwise,
(2.6)

where S(x, y; t) is the amount of signal S at location (x, y) at time t and α1
C and α2

C

are positive parameters.
The algorithms just described for the cell types in a BMU are summarized in

figure 4.

Figure 2.4: A proposal for BMU software. Signals R, S and T (in small circles)
are produced by the cell types listed. They may inhibit (-) or induce activation
(+) on the actions considered. Activation results from double inhibition, that is by
lowering the concentration of an inhibitor. Within any cell type, possible decision
choices are indicated by thin arrows. Dashed arrows correspond to a starting decision
stage in a newly formed cell. Releasing of signals by any cell type (for instance, R
and T in the case of active osteoblasts) is denoted by thick, white arrows.

Universidad Complutense de Madrid



2.3. MATHEMATICAL MODELLING 41

2.3.2 Spatial dependence in the model

In order to account for the effect of signals on cell decisions within a BMU, the spatial
disposition of both cells and molecules should be explicitly taken into account. To
do that, we will consider a 2 dimensional cross-section of bone adjacent to a section
of bone marrow by means of a cellular automaton (CA). More precisely, the bone
section considered will be thought of as a lattice divided in boxes of equal size, where
any such box can either remain empty or be occupied by only one individual cell.
In addition, boxes can contain different amounts of chemical cues. The dynamics of
signals in the CA is summarized as follows. At any given time, they are produced
by the corresponding cell types at constant rates. Signals are also assumed to un-
dergo Arrhenius-type decay (meaning that their decay rate is proportional to their
concentration) and to diffuse through the medium. We assume that diffusion occurs
faster that internal cell processes. In order to account for these two time scales si-
multaneously, we model diffusion by calculating, at each box and each iteration step
of the model, a weighted average of the amount of signals in neighboring positions
(see Figure 5.A). Cells need not remain confined to fixed boxes in the lattice. Active
osteoblasts and osteoclasts can decide to move anytime during the simulation. Such
decision is made upon comparing the amount of different signals in their adjacent
boxes; details about the movement of each cell type will be provided below. As a
starting point we consider a population of osteocytes (OCY), regularly distributed
within the bone matrix, and a layer of osteoblasts lining the bone section (see Figure
5.B). Simulations of the model are started after apoptosis of a group of osteocytes
has occurred (see Figure 5.C).

2.3.3 Cell movements

In this section we will describe the movements of active osteoblasts and osteoclasts
during a bone remodeling event.

Active osteoclasts

After osteocyte apoptosis, osteoclasts are recruited into the BRC by active os-
teoblasts. Once activated, osteoclasts star to remove old bone, thus giving raise
to the appearance of “cutting cones”, that may reach a variable depth into old bone
matrix, depending on the number of apoptotic osteocytes [22, 23]. We will assume
that osteoclasts remove bone at a rate that depends on the amount of signal S, thus
simulating the inhibitory effect that has been described for signals produced by os-
teocytes of osteoclasts activity [23]. We will assume that osteoclasts move to and
adjacent box once they have removed bone in their current position.

Bone resorption will be modeled by means of the following equations:b′C(x, y; t) = −αS3

C (1− S(x, y; t)

αS4

C

) if S(x, y; t) ≤ αS4

C

b′C(x, y; t) = 0 otherwise,

(2.7)

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 42

Figure 2.5: Spatial aspects of the model. A) Signal diffusion takes place at a
different scale than cell decisions. Let ∆t be the characteristic time step for cell
processes. In order to model signals diffusion, we calculate the profile resulting from
cells diffusing during this time interval. Diffusion can then be modeled at this char-
acteristic time by applying the corresponding profile to any source of signals in the
CA. B) Snapshots of a signal diffusion from the source shown in A. Time increases
from left to right and from top to bottom. C) Spatial arrangement of the main
elements of the model. Bone matrix is locally represented by means of a lattice con-
sisting of equal-size boxes, each able to accommodate one cell. The lower boundary
of the box considered in B is represented as lined with an osteoblast layer, lying
upon a layer of osteoclast precursors. D) Initial configuration of the model. At a
selected place within the bone matrix, osteocytes (blue dots) may undergo apoptosis,
thus triggering BMU operation in the corresponding region (blank). Lattice element
considered in B) is represented her for comparison purposes (see lower left corner).
Shades of blue represent different concentrations of signals in each box.

where S(x, y; t) is the amount of signal S at location (x, y) at time t and αS3

C and

αS4

C are positive parameters.

Active osteoblasts

Active osteoblasts move from their original position in the interface between bone
matrix and bone marrow to the inner surface of the cutting cone in order to secrete
new osteoid matrix. Osteoid production is also accompanied by active osteoblasts
movement, this time in the opposite direction. In our model, the movement of active
osteoblasts will be modeled as a consequence of cell division in the first stages of the
bone formation phase. The production of osteoid will, in turn, entail the apoptosis

Universidad Complutense de Madrid



2.4. RESULTS 43

and differentiation of some of the activated osteoblasts, which will result in the
movement of the active osteoblast layer back to the position occupied by osteoblasts
before the remodeling event (see Appendix A and figures A.1 and A.2 for details).

2.4 Results

In this Section we present some results obtained upon simulation of the model whose
elements have been described in our previous Section. To keep the flow of the main
arguments here, we postpone some technical details to Appendix A at the end of
the paper. Specifically, the reader will find there a flow chart describing the order in
which the corresponding algorithms are applied, as well as a list of the parameters
used in the simulations whose results will be shown below. Concerning this last issue,
we wish to emphasize that no parameter-fitting attempt has been made here. Our
concern is to show that the model proposed has not only the potential to reproduce
standard features in bone remodeling, but also to draw conclusions about the manner
in which this process is carried out that could not be a priori anticipated before
simulating that model. We will thus describe a choice of parameters for which such
goals are achieved, although we do not claim that these values are actual ones.
Parameters appearing in the actual BMU software should be of a structural nature,
and their precise values have to be determined experimentally. In fact, the circuitry
underlying BMU operation is expected to be rather sensitive to small variations
in such parameter values, a fact repeatedly observed in processes selected during
evolution. An example of such situation is provided by the Krebs cycle [42, 43]
which consists in a network of biochemical reactions taking place exactly at very
precise rates and proportions. We next succinctly describe some of the conclusions
drawn from the simulation of the algorithms described in Appendix A.

2.4.1 The BMU software proposed successfully recapitulates
bone remodeling

Upon induction of osteocyte apoptosis in a region, the cell algorithms proposed
sequentially describe the manner in which bone remodeling occurs there. See Figure
6 below.

2.4.2 The size of the BRC adapts to the number of apoptotic
OCY

The model proposed can be run when OCY apoptosis is assumed to occur in regions
of different sizes. In each case, the area of the corresponding Bone remodeling Com-
partment (BRC) is shown to depend on the extent of the induced OCY apoptosis.
See Figure 7.

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 44

Figure 2.6: Snapshots corresponding to sequential stages in a bone remod-
eling event. A ) A planar bone region is considered where active osteocytes (OCY,
blue dots) are regularly distributed. Osteoblasts (OBL, in green) are located at the
lower side, which represents the interface between bone and bone marrow. B) OCY
apoptosis is induced in a subset (in white) of the previous region, which leads to a
drop in OCY inhibitory action in that place. C) Lack of enough osteocyte inhibition
results in OBL activation, and the recruitment of osteoclast precursors which differ-
entiate into active osteoclasts, OCL (in red). D) OCL form cutting cones that move
into the apoptotic region destroying old bone (bone resorption). E) OCL- driven
resorption proceeds until the inhibitory effect of the remaining osteocytes precludes
further OCL progression. F) Activated osteoblasts (OBA, in orange) move in the
wake of OCL. G) OBA move inward from the OCL- lined cutting cone boundary
secreting osteoid matrix. Part of them subsequently differentiate into OCY (blue
dots). The limit of the cutting cone is shown in black. H) When remodeling is
finished, a new distribution of osteocytes is achieved in the remodelled region, and
I ) new bone is eventually formed.

2.4.3 Osteoclast lifespan depends on the depth of the BRC

A prediction of the model is that average osteoclast lifespan is not a priori deter-
mined, but depends instead on the size of the region where bone remodeling has to
be performed. This fact is illustrated in Figure 8 below, where mean OCL lifespan
is plotted in terms of BRC depth.

2.4.4 The BMU software is robust to signal noise

As we have seen above, the operation of cells within a BMU depends on cell-to-
cell communication by means of signals that diffuse through the bone matrix. The
coupling of old bone resorption and new bone formation can therefore be affected by
signal noise due, for instance, to bone matrix spatial heterogeneity. In figure 9 we

Universidad Complutense de Madrid



2.4. RESULTS 45

Figure 2.7: Simulation of three different scenarios of OCY apoptosis. Im-
ages correspond to results obtained in each case after the simulation of the model
since OCY apoptosis until the end of the remodeling process. The edge of the re-
sulting cutting cone is depicted in black. The shape and size of the resulting BRC
is determined in each case by the region where OCY apoptosis has occurred.

Figure 2.8: Dependence of OCL lifespan on BRC size. Left) A region where
bone resorption has occurred following OCY apoptosis is depicted. BCR depth is
defined there as the longest distance within that region measured perpendicular to
its lower side. Right) A plot of average OCL lifespan during bone remodeling events
in terms of BRC depth obtained after simulations of the model (each dot corresponds
to a simulation). Note that OCL live longer when BRC is larger.

show that the suggested BMU software is robust to environmental noise.

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 46

Figure 2.9: Effect of signal noise in BMU operation. A) Simulation cor-
responding to a bone remodeling without noise in signal S diffusion. The initial
situation is depicted in the left figure. The shape of the cutting cone is marked in
black. The region of apoptotic OCYs is shown in the center. The right picture cor-
responds to the new bone after remodeling. B) In spite of noise, both the resulting
cutting cone and the new distribution of OCY is similar to A. C) Histogram of the
rate between the number of newly formed OCY in the BRC after bone remodeling
with and without noise (the result was obtained after 200 simulation of the model).

2.5 Discussion

In this work we have proposed and analyzed a minimal model for bone remodeling
carried out by Bone Multicellular Units (BMUs) as part of their routine activity
to replace old bone by new one, a process which is often triggered by changes in
mechanical load. We have already recalled that such remodeling is continuously
taking place during human life, so that on average the whole skeleton is renewed
every 10 to 15 years [44]. BMUs activity is usually carried out at small regions, say
about 4000 microns in width [45] and therefore the total number of cells involved
in each remodeling event is comparatively small. The previous length scale can be
also used to characterize micro-fractures (arising for instance from physical exercise)
which are also self-repaired in the same manner. BMUs operation can thus be viewed
as the basic building block to keep homeostasis in bone tissue. When normal function
in bone is compromised by large-scale hazards (say, fractures), more sophisticated
repair mechanisms, involving in particular inflammatory signals, blood clot formation
and the production of different callus templates for bone repair are involved, which
have not been discussed here.

Our goal in this work is to identify, by means of a mathematical model, a minimal
algorithm for BMU operation. To do this, we have kept to a few basic principles.

Universidad Complutense de Madrid



2.5. DISCUSSION 47

First, at any time during that process, each cell within the different lineages present
in a BMU should individually decide one among a restricted set of actions: divide,
die or differentiate. Such choice is not assumed to be random, but determined
by feedback from the surrounding medium into the internal dynamics of molecules
which act as decision inhibitors. No previous planning is presupposed, genetically
or otherwise. Second, only a minimal number of signals are selected that allow for
BMU functioning. As recalled in Section 2, all elements retained have been shown
to be actually in place in working BMUs. In particular, cell decision inhibitors have
been identified, and effects induced by the signals proposed have been shown to be
caused by chemical mediators experimentally detected. However, we have not aimed
at producing a comprehensive model where all data currently known could be fitted.
We have instead attempted to show what the basic ingredients of an operational
BMU could be like. Implicit in this approach is the assumption that some signals
already identified are multi-functional (thus yielding different effects in different cell
lines ) and that in general some signaling networks observed can be redundant,
perhaps as a consequence of having been arrived at through different paths in the
course of evolution. Notice, however that we do not claim that our proposal is the
only possible one. In principle, there could be alternative minimal models for what
we have called a BMU software. We believe however, that the degree of complexity
retained in our model could not be significantly reduced by alternative mechanisms.

Selecting a minimal model has the advantage of dealing with a comparatively
small number of parameters. Indeed, it is well known that a given phenomenon can
be fitted by utterly different models, relying on different principles, upon appropriate
choices of a suitable (and often rather modest) number of parameters appearing in
them [46, 47]. We have shown herein that the parameters retained in our model
can be selected so that BMU standard operation is reproduced and, importantly,
consequences such as the dependence of osteoclasts lifespan on the size of the region
to be repaired can be inferred that were not a priori clear before the model was sim-
ulated. We do not claim, however, that the results obtained provide a validation of
the model. In our opinion, such validation could only stem from experimental deter-
mination of such values. This should be in principle possible, since most parameters
retained are of a structural character (see pseudo-code guidelines in Appendix A ).
By that we mean that such values correspond to coefficients as for instance exter-
nal signals/internal inhibitors interaction or maximal inhibitor accumulation within
osteoblasts. Such performance rates might well have been evolutionary selected to
remain confined within very narrow margins, as is the case with virtually all robust
metabolic processes which have been unraveled so far [42, 43].

We conclude this discussion by briefly remarking on possible future directions
arising from this work. To begin with, the experimental question of determining
the operational parameters introduced in our model should be again stressed. If
and when this is done, we would have in place a BMU chip structure to be used as
a building block for operational purposes. On the other hand, one could simulate
on models as that presented here the impact of external cues in the modulation
of the repair process, either to slow it down or to accelerate it. Finally, as the

Universidad Complutense de Madrid



CHAPTER 2. BONE REMODELLING 48

size of the region to be repaired increases, a threshold should be reached beyond
which inflammatory signals become relevant, and a new and more complex stage is
entered. Understanding the matching between these two cases could be instrumental
in designing techniques to foster bone fracture repair. We intend to pursue some of
these subjects elsewhere.

Universidad Complutense de Madrid



Bibliography

[1] J. M. Graham, B. P. Ayati, S. A. Holstein and J. A. Martin (2013). The role
of osteocytes in targeted bone remodeling: a mathematical model. PLoS ONE
8(5): e63884. doi:10.1371/journal.pone.0063884

[2] M. L. K. Tate, J. R. Adamson, A. E. Tami and T. W. Bauer (2004). The
osteocyte. Int. J. Biochem. Cell Biol. 36, 1: 1- 8.

[3] F. F. Safadi, M. F. Barbe, S. M. Abdelmagid, M. C. Rico et al. (2009). Bone
structure, development and bone biology. In Bone Pathology; 1.50. Humana
Press.

[4] R. L. Jilka, R. S. Weinstein, A. M. Parfitt and S. C. Manolagas (2007). Perspec-
tive : Quantifying osteoblast and osteocyte apoptosis: challenges and rewards.
Journal of Bone and Mineral Research 22, 10: 1492-1501.

[5] B. S. Noble, N. Peet, H.Y. Stevens, A. Brabbs et al. (2003). Mechanical loading:
biphasic osteocyte survival and targeting of osteoclasts for bone destruction in
rat cortical bone. Am. J. Phys-Cell Phys. 284, 4: C934-C943.

[6] J. P. Bidwell, J. Yang and A. G. Robling (2008). Is HMGB1 an osteocyte alarmin
?. J. Cell. Biochem. 103, 6: 1671-1680.

[7] A. M. Parfitt (2002). Targeted and nontargeted bone remodeling: relation to
basic multicellular unit and progresi?n. Bone 30(1): 5-7.

[8] E. F. Eriksen (2010). Cellular mechanisms of bone remodeling. Reviews in en-
docrine and metabolic disorders 11(4): 219-227.

[9] E. M. Hauge, D. Qvesel, E. F. Eriksen, L. Mosekilde and F. Melsen (2001).
Cancellous bone remodeling occurs in specialized compartments lined by cells
expressing osteoblastic markers. J. Bone. Min. Res.16, 9: 1575- 1582.

[10] I. H. Kalfas (2001). Principles of bone healing. Neurosurgical Focus 10, 4: 1-4.

49



BIBLIOGRAPHY 50

[11] T. A. Franz-Odendaal, B. K. Hall and P. E. Witten (2006). Buried alive: how
osteoblasts become osteocytes. Dev. Dyn. 235, 1: 176-190.

[12] N. A. Sims and T. J. Martin (2014). Coupling the activities of bone forma-
tion and resorption: a multitude of signals within the basic multicellular unit.
BoneKEY reports 3.

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

[14] M. Ryser, S. Komarova and N. Nigam (2010). The cellular domains of bone
remodeling: a mathematical model. SIAM J. Appl. Math. 70: 1899-1921.

[15] J. M. Graham, B. P. Ayati, S. A. Holstein and J. A. Martin (2013). The role
of osteocytes in targeted bone remodeling: a mathematical model. PLoS ONE
8(5), e63884.

[16] A. G. Robling, A. B. Castillo and C. H. Turner (2006). Biomechanical and
molecular regulation of bone remodeling. Annu. Rev. Biomed. Eng. 8: 455-498.

[17] K. T. Steeve, P. Marc, T. Sandrine, H. Domenique and F. Yannick (2004).
IL-6, RANKL, TNF-α/IL-1: interrelations in bone resorption pathophysiology.
Cytokine & growth factor reviews 15 (1): 49-60.

[18] B. S.Noble (2008).The osteocyte lineage. Archives of Biochemistry and Bio-
physics 473, 2: 106-111.

[19] G. Chen, C. Deng and Y. P. Li (2012). TGF- β and BMP signalling in osteoblast
differentiation and bone formation. Int. Journal of Biological Sciences 8, 2: 272.

[20] R. Marsell and T. A. Einhorn(2011). The biology of fracture healing. Injury 42,
6: 551-555.

[21] X. Cao (2011). Targeting osteoclast- osteoblast communication. Nature
Medicine 17,11: 1344.

[22] T. J. Heino, T. A. Hentunen and H. K. Väänänen (2002). Osteocytes inhibit
osteoclastic bone resorption through transforming growth factor- β? Enhance-
ment by estrogen. J. Cell. Biochemistry 85, 1: 185-197.

[23] G. S. Stein, J. B. Lian and T. A. Owen (1990). Relationship of cell growth to
the regulation of tissue-specific gene expression during osteoblast differentiation.
The FASEB Journal 4, 13: 3111-3123.

[24] J. M. Hock, V. Krishnan, J. E. Onyia, J. P. Bidwell, J. Milas and D. Stanislaus
(2001). Osteoblast apoptosis and bone turnover. Journal of Bone and Mineral
Research 16, 6: 975-984.

Universidad Complutense de Madrid



BIBLIOGRAPHY 51

[25] W. J. Boyle, W. S. Simonet and D. L. Lacey (2003). Osteoclast differentiation
and activation. Nature 423 (6937): 337-342.

[26] S. L. Teitelbaum (2000). Bone resorption by osteoclasts. Science 289 (5854):
1504-1508.

[27] B. S. Noble and J. Reeve (2000). Osteocyte function, osteocyte death and bone
fracture resistance. Molecular and Cellular Endocrinology 159, 1: 7-13.

[28] C. A. Mullen, M. G. Haugh, M. B. Schaffer, R. J. Majeska and L. M. McNamara
(2013). Osteocyte differentiation is regulated by extracellular matrix stiffness
and intercellular separation. J. Mech, Behav. Biomed. Mat. 28: 183-194.

[29] E. Abe (2006). Function of BMPs and BMP antagonists in adult bone. Ann.
New York Acad. Sci. 1068, 1: 41 - 53.

[30] M. J. Christopher and D. C. Link (2008). Granulocyte Colony-Stimulating Fac-
tor induces osteoblast apoptosis and inhibits osteoblast differentiation. Journal
of Bone and Mineral Research 23, 11: 1765-1774.

[31] L. D. You, S. Weinbaum, S. C. Cowin and M. B. Schaffer (2004). Ultrastructure
of the osteocyte process and its pericellular matrix. The Anatomical Record A.
278, 2: 505-513.

[32] J. P. Reyes, S. M. Sims and S. J. Dixon (2011). P2 receptor expression, signalling
and function in osteoclasts. Front. Biosci. (Schol. Ed) 3: 1101-1118.

[33] M. Karin and T. Hunter (1995). Transcriptional control by protein phospho-
rylation: signal transmission from the cell surface to the nucleus. Curr Biol
5:747757.

[34] R. L. Jilka, R. S. Weinstein, T. Bellido, A. M. Parfitt and S. C. Manolagas
(1998). Osteoblast programmed cell death (apoptosis): modulation by growth
factors and cytokines. Journal of Bone and Mineral Research, 13(5), 793-802.

[35] C. F. Arias, M. A. Herrero, F. J. Acosta and C. Fernández -Arias (2014). A
mathematical model for a T cell fate decision algorithm during immune re-
sponse. Journal of Theoretical Biology 349: 109-120.

[36] P. Ducy (2000). Cbfa-1: a molecular switch in osteoblast biology. Dev. Dyn.
219,4 ; 461-471.

[37] P. Bialek, B. Kern, X. Yang, M. Schrock et al. (2004). A twist code determines
the onset of osteoblast differentiation. Dev. Cell 6, 3 ; 423-435.

[38] M. Zhao, M. Qiao, B. Oyajobi, G. R. Mundy and D. Chen (2003). E3 ubiquitin
ligase Smurf 1 mediates core-binding factor ?1/Runx2 degradation and plays a
specific role in osteoblast differentiation. J. Biol. Chem. 278, 30: 27939-27944.

Universidad Complutense de Madrid



BIBLIOGRAPHY 52

[39] W. Huang, S. Yang, J. Shao and Y.P. Li (2007). Signaling and transcrip-
tional regulation in osteoblast commitment and differentiation. Frontiers in
Bioscience: 12 ,

[40] H. M. Ryoo, M. H. Lee and Y. J. Kim (2006). Critical molecular switches
involved in BMP-2 induced osteogenic differentiation of mesenchymal cells.Gene
366, 1: 51-57.

[41] J.F. Chau, W. F. Leong and B. Li (2009). Signaling pathways governing os-
teoblast proliferation, differentiation and function. Histology and histopathology,
24(12): 1593-1606.

[42] D. L. Nelson, A. L. Lehninger and M. M. Cox.Lehninger principles of biochem-
istry. Macmillan, 2008.

[43] O. Ebenhöh and R. Heinrich (2001). Evolutionary optimization of metabolic
pathways. Theoretical reconstruction of the stoichiometry of ATP and NADH
producing systems. Bulletin of mathematical biology, 63(1): 21-55.

[44] M. Zaidi (2007). Skeletal remodeling in health and disease. Nat. Med. 13: 791-
801.

[45] A. M. Parfitt (2002). Targeted and non targeted bone remodeling: relationship
to basic multicellular unit origination and progression. Bone 30, 1: 5-7.

[46] H. C. Hemker, S. Kerdelo and R. M. W. Kremers (2012). Is there value in kinetic
modeling of thrombin generation? J. Thromb. Haemost. 10: 1470-1477.

[47] L. R. Ginzburg and C. X. J. Jensen (2004) . Rules of thumb for judging ecological
theories. Trends in Ecology and Evolution 19, 3 ; 121-126.

Universidad Complutense de Madrid



Chapter 3
Early stages of Bone fracture healing

3.1 Introduction

Bone repair is a robust physiological process which is constantly acting in healthy
individuals. In particular, bone microfractures induced by exercise (including walk-
ing) are routinely self-healed without even being noticed, although large microfrac-
ture condensation may result in severe injury [39], [6]. Interestingly, bone self-repair
makes use of key ingredients of the developmental program responsible for ossifica-
tion during the embryonic stage [19], [18]. The sequence of events occurring after a
bone fracture involves three phases: inflammatory, reparative and remodelling, each
characterized by specific cellular, histological and biochemical features. The inflam-
matory phase has received relatively lower attention then the other two, in spite of
its fundamental role in the subsequent evolution of the process.

In this chapter we propose and discuss a simple model to analyze the early stages
of bone fracture repair. When a bone fractures, several blood vessels break and a
mass of clotted blood is formed at the injury site. This fracture hematoma has a
little effect to stabilize the bone, but it provides a first bridge between the two bone
fragments and thus keeps both pieces lined up for mending. A key ingredient in this
stage is the formation of the first elastic, fibrin-collagen scaffold (the so-called early
fibrous callus, EFC) which provides a first consistent scaffold for mesenchymal cells
to settle and differentiate to fibroblasts, condrocytes and osteoblasts, to eventually
form a bony bridge between the two fragments. The corresponding process, which is
an essential starting point for subsequent bone healing, takes place in a comparatively
short period: no more than 48 hours in the case of thin induced fractures in mice [21].
More precisely, we shall pursue the following goals: (i) to estimate the characteristic
times of the main events leading to EFC formation, and (ii) to ascertain their
dependence on the kinetics of key underlying mechanisms as platelet activation, fibrin
polymerization and fibrin clot remodelling mediated by fibroblasts. We will not be
directly concerned here with mechanical effects, which become particularly relevant

53



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 54

at later stages, as EFC becomes a cartilaginous callus (CC) on which different types
of bone tissues are formed upon remodelling of CC through a sequence of processes
which result in changes in the composition and mechanical properties of previous
scaffolds [21], [20].

Our plan is as follows. In Section 2 below we review basic known biological facts
concerning the main events in fracture repair, particularly during the early stage to
be considered. We then discuss in Section 3 a mathematical model whose analysis
will allow us to estimate several characteristic times of the process under consider-
ation . The first of them corresponds to the formation of a fibrin blood clot, the
earliest structure generated during the process of bone repair. The onset of such a
clot, which in spite of its poor mechanical properties provides a first bridge to join
both sides of the fracture and is therefore crucial for bone repair [50] is character-
ized here as a phase transition in a polymerization event. Its timing is estimated
in terms of the various physiological processes (thrombin-induced activation, fibrin
polymerization,...) involved.

Once this fibrin network is in place, attention is paid to a first , and crucial,
remodelling effect. This results from the release of collagen over the fibrin net,
which is deposited by fibroblasts arising from asymmetric division of mesenchymal
stem cells (MSCs) recruited from neighbouring regions (periostium). The timing of
this event is kept track of by means of the motion of a front of MSCs proceeding
inwards from the borders of the fracture gap. Such motion is chemotactically induced
by growth factors released by activated platelets trapped within the fibrin clot [54],
[24]. As it turns out, substantial insight about the process considered is derived from
the study of a particular asymptotic limit, namely the case where fibrin clotting is
triggered by an instantaneous thrombin discharge. Under such assumption, our
model can be explicitly solved, which permits a detailed discussion of the influence
of the different biological processes involved. Estimates for the characteristic times
obtained in Section 3 are then provided in Section 4. Finally, a discussion on the
results obtained as well on possible future directions is gathered in a concluding
Section 5.

We point out that mathematical models of different aspects of ossification and
bone repair have received considerable attention in recent years [35],[2],[23],[16],[44];
for later repair stages see also [41], [40], [45]. Some issues similar to those addressed
in this work have been also discussed in related biological problems as wound healing,
which shares a number of features with bone repair [38], [42], [32] in particular during
the early stages here considered.

3.2 Basic facts on bone fracture healing. The first
48 hours

When bone fracture occurs, damage is induced on a variety of tissues including
surrounding muscular cells and periosteum as well as on bone itself . Bone is a
highly vascularized tissue, and its fracture leads to many blood vessels acrossthe

Universidad Complutense de Madrid



3.2. BASIC FACTS ON BONE FRACTURE HEALING. THE FIRST 48 HOURS 55

fracture line being broken . Blood then runs away from disrupted vessels to the
medullary canal between the fracture ends and rapidly coagulates to form a clot. On
the other hand, immediately after fracture the resulting necrotic material induces an
intense inflammatory response including widespread vasodilation, plasma exudation,
edema and recruitment of acute inflammatory cells [9]. As a consequence, a complex
hematoma is generated that is of major importance for the beginning of bone healing,
since it provides a first bridge between separated fragments of bone [50]. Such
hematoma plays a very small mechanical role in immobilizing the fracture, but it
serves as a fibrin scaffold and provides an adequate microenvironment over which
repair cells may settle to perform their function.

Formation of hematoma is a comparatively fast process. In the case of thin in-
duced fractures in rat tibia, a case which will be considered here for comparison
purposes, the fracture gap is about 2.5 mm in average, and hematoma can be ob-
served by optical microscopy at about 3 hours after fracture induction. See Figure
1 below.

Figure 3.1: Blood clot between the ends of a bone fracture. Bone tissue is stained in pink
whereas blood tissue has a brownish color. The immediate ends of the fracture border contain
necrotic material (arrows). Picture take from tibia from a 35 days old rat, 6 hours after fracture.

A key constituent of hematoma is a fibrin net. Immediately after bone breaking,
a coagulation process starts. Thrombin generated by blood vessel disruption induces
the cleavage of fibrinogen present in blood to yield fibrin monomers, which quickly
polymerize afterwards. In this manner a fibrin net is generated that in 6 to 8 hours
connects both sides of the fracture and provides a first scaffold on which bone will
eventually be produced [21], [15]. Importantly, when fracture occurs, a number of
growth factors ( BMP, VEGF, PDGF, .. ) are released within the fracture gap
[5]. Such substances will mediate subsequent steps in fracture repair.

As soon as a fibrin clot begins to form, new functional blood vessels as well as spe-
cific cell populations are recruited from the fracture boundaries. The cells involved
directly in the repair of fractures are of mesenchymal origin and of a pluripotential
nature. For definiteness we shall refer to such precursors as mesenchymal stem cells
(MSCs) in the sequel. MSCs move in the wake of the new vascular bed created by
incoming blood vessels. The latter can be seen within 6 to 12 hours near the fracture
boundaries and between 24 to 48 hours they have invaded the whole fracture gap

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 56

[51]. See Figure 2 below where images taken 24 hours after fracture induction are
presented.

Figure 3.2: Early remodelling can be perceived in hematoma 24 hours after fracture, with fibrillar
material being found in some regions of the blood clot represented in pink (arrows). Bone depicted
in brown.

Figure 3.3: Image showing early invading cells derived from the cambium layer of the periosteum
at the border of the fracture 24 hours after fracture. Such cells show morphology of mesenchymal
cells (arrows).

Figure 3.4: Picture representing a EFC being formed 24 hours after fracture. An increasing
number of invading cells is observed (arrows) and partial reabsorption of the hematoma can be
noticed (stars)

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 57

MSCs are chemotactically recruited from the fracture borders by chemical cues
released from activated platelets trapped in the fibrin clot [30] (see Figure 3.3).
MSCs may sustain either symmetric or asymmetric division [36]. Symmetric di-
vision, whereby a MSC gives raise to two new MSCs, takes place whenever their
initial number is not enough to keep the healing process in time. On the other hand,
when a MSC reaches the fibrin net, it can undergo asymmetric division, thus giving
raise to one MSC and one differentiated cell, which may be of various types: fibrob-
lasts, chondrocytes or osteoblasts depending on microenvironment stimuli and on
the stresses they are subject to [36] [61], [31]. In particular, such lineages are known
to depend on the density of the net on which division takes place [11], [14], [48]. Cell
replication is stimulated by some of the growth factors present, whose action is in
turn dependent on the unfolding of the fibrin clot. In particular, anchorage into the
net significantly extends the half-life of growth factors with respect to that of freely
diffusing ones [49], [60].

Specifically, the first type of asymmetric division takes place over a low-density
fibrin clot and gives raise to fibroblasts. Once attached to the net, fibroblasts secrete
a number of substances (collagens) which stick to the clot and change its mechanical
properties [14], [48]. The resulting network, called early fibrous callus (EFC) is in
place between 24 to 48 hours after bone break up (see Figure 3.4). It has larger
mechanical stability than its preceding fibrin net, although not enough to ensure
working stability to the healing bone. Such stability will be achieved through subse-
quent remodeling processes. In particular, upon sticking to a collagen-enriched fibrin
net, MSC can again undergo asymmetric division that now gives raise to one MSC
and one chondrocyte [36], [11]. This new type of cell interacts with the existing net
and produces a more stable compound, called cartilaginous callus. Bone formation
properly begins in a subsequent remodeling phase. At that stage, cartilaginous cal-
lus is replaced first by trabecular and then by lamellar bone [55] through a process
similar to endochondral ossification. These later stages will not be considered in this
work.

3.3 A mathematical model for the formation of a
fibrin-collagen callus template.

In this Section we formulate and analyze a mathematical model to describe the early
stages of cartilaginous callus formation in a fracture. For simplicity we specialize
here to a one-dimensional case, and assume that the fracture gap is located at an
interval centered around the origin, say Ω =

(
−L

2 ,
L
2

)
, with L > 0. Our model is

of a modular character, and contains a series of basic units which are activated in
sequential order. These will be described below.

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 58

3.3.1 Generation of fibrin monomers upon thrombin-induced
fibrinogen cleavage.

Right after a fracture is generated, a number of substances arising from broken
blood vessels are allowed to enter in contact with leaked blood. The actual sequence
of biochemical events thus triggered is known to involve several agents [17], [43].
However, a common simplifying assumption consists in considering that complex
process as driven by the action of the enzyme thrombin (obtained from its inactive
precursor prothrombin) which cuts inactive fibrinogen molecules present in blood,
thus giving raise to fibrin monomers. The latter then quickly polymerize to produce
a fibrin clot. It is customary to describe thrombin generation by means of curves of
the form represented in Figure 3.5, below ([17], [8]).

Thrombin 

Concentration

Time (Min)

(A)

Thrombin 

Concentration

Time (min)(B)

Figure 3.5: Typical thrombin profiles (adapted from [17], [8]). (A) About 95% of thrombin is
generated as a consequence of positive feedback processes involving coagulation factors FVIII,
FIX and FX [17], [37]. (B) dependence of thrombin concentration on time in some experimental
settings, when anticoagulant factors are kept below haemostatic levels [8].

As mentioned in the caption of Figure 3.5, the sequence of events leading to
thrombin generation include a number of feedback loops of an activator-inhibitor
nature involving several coagulation/anticoagulation factors ([17], [25], [27]). Given
that large uncertainties remain about the order and rates of the corresponding chem-
ical reactions ([12], [57]), we have resorted here to a simplified modelling approach,
which is now described. Let us respectively denote by F (t), θ(t) the concentration of
fibrinogen and thrombin (assumed homogeneous in space) at time t. We postulate
that fibrinogen is removed by thrombin according to a second order kinetics

Ḟ = −kgFθ for some kg > 0, t > 0 (3.1)

whereas at t = 0, when fracture occurs, one has that

F (0) = FM > 0 (3.2)

for some physiological concentration FM . We now consider an asymptotic limit of
the profile sketched in Figure 3.5 A, namely the case of an instantaneous thrombin
discharge at t = TM > 0, that is:

θ(t) = Cδ(t− TM ) for some C > 0 (3.3)

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 59

where δ(t − TM ) denotes Dirac’s delta at t = TM . Note that (3.1) and (3.3) rep-
resent a very fast chemical reaction for which the precise mathematical meaning of
the governing differential equation has to be specified. In particular (3.1) is to be
understood here as stating:

d

dt
(logF ) = −kgθ

On integrating this equation and using (3.2) it follows at once that:

F (t) =

{
FM for t < TM

FMe−kgC for t > TM
(3.4)

Fibrinogen molecules cleaved by thrombin give rise to fibrin monomers, which
subsequently polymerize to form large fibrin agregates. Using classical polymeriza-
tion theory ([53], [28], [25]) the corresponding process can be modelled by an infinite
system of coagulation-fragmentation equations for the concentrations Fk(t) of fibrin
k-mers k ≥ 1

Ḟ1 = kgFθ − F1

∞∑
i=1

a1iFi +
∞∑
i=1

b1iFi+1 − krF1 (3.5)

Ḟk =
1

2

k−1∑
m=1

am,k−mFmFk−m − Fk

∞∑
i=1

akiFi +
∞∑
i=1

bkiFi+1

− 1

2

∑
i+j=k

bijFk − krFk; k = 2, 3, . . .

(3.6)

where akj , (respectively bkj) denote the coagulation rates between k-mers and j-mers
(respectively the fragmentation rate at which k-mers are produced from the break-up
of j-mers, j > k). Rates kr correspond to polymer deactivation, and will be assumed
to be independent of polymer length by simplicity. Assuming polymerization to
occur according to classical Flory-Stockmayer theory, so that no polymer rings can
be formed and coagulation rates are proportional to the number of the free available
sites in each polymer, it turns out that

aij = kp ((f − 2)i+ 2) ((f − 2)j + 2)

([25], [28]) where f denotes the number of active functional sites in a fibrin monomer.
Taking f = 4 ([59]) we eventually obtain

aij = 4kp(i+ 1)(j + 1) (3.7)

On the other hand, assuming (i) the probability of a molecule of length k being cut
by a fragmentation mechanism to be proportional to the number (k − 1) of bonds
and (ii) the probability to form a molecule of length i out of fragmentation of a
molecule of length k (k > i) to be inversely proportional to the number of bonds, we
obtain ([25]) that:

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 60

bij = 2kb for some kb > 0 (3.8)

A general discussion of the resulting system (3.5)-(3.8) is highly involved and will
not be pursued here. However, considerable insight into blood-clotting dynamics is
gained by making use of suitable auxiliary assumptions and simplifying hypotheses.
To this end, for any integer n ≥ 0 we define the n-moment Mn of the polymer
distribution:

Mn(t) =

∞∑
k=1

knFk(t) (3.9)

Note that M1(t) represents the total mass of the polymer distribution once the
monomer mass has beeen normalized to one. Of particular interest here is the fact
that the ratio (M2/M1) corresponds to the average molecular size. In classical Flory-
Stockmayer theory a finite-time blow-up of that ratio marks the onset of a sol-gel
phase transition, whereby a gel phase appears ([53]). Bearing this fact in mind,
we say that a fibrin clot is formed as soon as that ratio develops a singularity at
some time ([28], [25]). Upon replacing kb by 3kb, it follows from (3.5)-(3.8) that the
moments M1, M2 satisfy:

Ṁ1 = kgFθ − krM1

Ṁ2 = kgFθ + 4kp (M2 +M1)
2 − kb(M3 −M1)− krM2

(3.10)

Note that such a system is not in a closed form, since the third moment M3 appears
in equation for M2. To obtain a well posed problem we now set:

M3 =
M2

2

M1
(3.11)

As observed in [25] where this closure argument was proposed, (3.10) holds true in
two extreme cases, namely when either only monomers exist or when all fractions are
concentrated in a single, large cluster and no smaller fractions are present. However,
one always has that M2

2 ≤ M1M3 by Holder’s inequality. Thus, when substituting
(3.11) into (3.10), equality sign in the second equation for (3.10) has to be replaced
by ≤ there. Therefore, the closed system

Ṁ1 = kgFθ − krM1 (3.12)

Ṁ2 = kgFθ + 4kp (M2 +M1)
2 − kb

(
M2

2

M1
−M1

)
− krM2 (3.13)

is satisfied by solutions M1 to (3.12) and by subsolutions M2 to (3.13).

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 61

3.3.2 Estimating fibrin clot formation

As it turns out, equations (3.12), and (3.13) can be explicitly solved under our current
assumption (3.3). Notice that using (3.1), (3.12) can be rewritten in the form

d

dt

(
ekrtM1

)
= kge

krtFθ = −ekrtḞ

it then follows that

M1(t) = 0 for t < TM ; M1(TM ) = FM

(
1− e−kgC

)
≡ A (3.14)

M1(t) = Ae−kr(t−TM ) for t > TM (3.15)

Concerning the second moment, one also has that

M2(t) = 0 for t < TM ; M2(TM ) = A (3.16)

Notice that the values M1(TM ) = M2(TM ) = A are obtained from the fact that by
(3.1), (3.2) and (3.13) the jump of M1, M2 at t = TM is provided by instantaneous
activation of fibrinogen by thrombin at that time. We next point out that, on
multiplying both sides in (3.13) by e−kr(t−TM ) and upon noticing that e−kr(t−TM ) =
M1(t)

A for t > TM (cf. (3.15)) one readily sees that U(t) = M2(t)
M1(t)

satisfies:

dU

dt
= 4kpAe−kr(t−TM )(U + 1)2 − kb(U

2 − 1) for t > TM (3.17)

U(TM ) = 1. (3.18)

Equation (3.17) is of Riccati type. Since U(t) = −1 is a particular solution, its
general solution can be sougth in the form

U(t) = Y (t)− 1 (3.19)

To better understand the consequences of the various kinetic processes retained in
(3.17), the following situations will be separately considered:

1. Case: kp > 0, kb = kr = 0. When fragmentation and deactivation are
considered negligible with respect to coagulation, we obtain from (3.17) and
(3.19) that

U(t) =
1 + 8Akp(t− TM )

1− 8Akp(t− TM )
(3.20)

which blows up at a time Tf given by

Tf = TM +
1

8Akp
(3.21)

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 62

2. Case: kp > 0, kb > 0, kr = 0.

We now obtain from (3.17) and (3.19) that

U(t) =
kb + 4Akp

(
1− e−2kb(t−TM )

)
kb − 4Akp

(
1− e−2kb(t−TM )

) (3.22)

U(t) blows up at t = Tf if

1− e−2kb(Tf−TM ) =
kb

4Akp
(3.23)

which requires
kb < 4Akp (3.24)

Therefore, large fragmentation rates preclude the formation of fibrin clots.
When (3.23) holds, it follows from (3.23) that

Tf = TM − 1

2kb
log

(
1− kb

4Akp

)
(3.25)

Note that (3.25) reduces to (3.21) in the limit kb → 0. Moreover Tf (kb)
is increasing in kb, so that (3.21) provides a lower bound for (3.25) when
0 < kb ≪ 1.

3. Case: kp > 0, kb > 0, kr > 0.

When all kinetic rates are different from zero substitution of (3.19) into (3.17)
leads to

dY

dt
= f(t)Y 2 + 2kbY ; f(t) = 4kpAe

−kr(t−TM ) − kb (3.26)

so that Z = − 1
Y solves

dZ

dt
+ 2kbZ = f(t), Z(TM ) = −1

2

whence

Z(t) =e−2kb(t−TM )

(
−1

2
+

∫ t

TM

e2kb(s−TM )f(s)ds

)
=

4Akp
2kb − kr

(
e−kr(t−TM ) − e−2kb(t−TM )

)
− 1

2

(3.27)

for t > TM when 2kb ̸= kr, and:

Z(t) = 4Akp(t− TM )e−2kb(t−TM ) − 1

2
(3.28)

for t > TM and 2kb = kr. Assume now that 2kb > kr then, on setting:

Q =
4Akp

2kb − kr
=

4FM (1− e−kgC)kp
2kb − kr

(3.29)

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 63

it follows that U(t) = M2(t)
M1(t)

is given by:

U(t) =
1 +Q

(
e−kr(t−TM ) − e−2kb(t−TM )

)
1− 2Q

(
e−kr(t−TM ) − e−2kb(t−TM )

) (3.30)

Let g(t) = e−kr(t−TM ) − e−2kb(t−TM ) for t ≥ TM . Clearly g(TM ) = 0, g(t) > 0
for t > TM and g(t) achieves a maximum at some time t = t∗. A straightfor-
ward computation yields that:

t∗ = TM +
1

2kb − kr
log

(
2kb
kr

)
(3.31)

m∗ = g(t∗) =

(
2kb
kr

)− kr
2kb−kr

−
(
2kb
kr

)− 2kb
2kb−kr

(3.32)

Thus for U(t) in (3.30) to blow up at time Tf , we need

2Qm∗ ≥ 1 (3.33)

where Q, m∗ are respectively given in (3.29) and (3.32). If (3.33) holds, one
has that

Tf ≤ TM +
1

2kb − kr
log

(
2kb
kr

)
(3.34)

the case where 2kb < kr is similarly dealt with. Finally, when 2kb = kr blow
up for U(t) occurs if there exists Tf > TM such that

8Akp(Tf − TM )e−2kp(Tf−TM ) = 1 (3.35)

It is readily seen that (3.35) holds if m∗ = maxt>TM

(
(t− TM )e−2kb(T−TM )

)
is such that 8Akpm

∗ ≥ 1. Since m∗ = (2ekb)
−1 in this case, the resulting

condition reads:

4Akp > ekb (3.36)

We conclude this section by observing that our discussion can be extended to the
case where (3.3) is replaced by

θ(t) = C(x)δ(t− TM )

for some function C(x). One just has to replace C by C(x) in the corresponding
estimates, as for instance in (3.21). Note that this implies that Tf = Tf (x) so that
different clot nucleation centers may be generated at different times.

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 64

3.3.3 Platelet activation and growth factor release

Besides cleaving fibrinogen, thrombin also mediates platelet activation [17]. Acti-
vation leads to platelets changing shape, sticking to each other and releasing the
content of their storage granules [7]. In the sequel we will denote by pa(t) the con-
centration of activated platelets at time t > 0. Ignoring deactivation effects, we
propose the dynamics of pa(t) to be given by:

ṗa(t) = α(P0 − pa(t))θ for t > 0 (3.37)

pa(0) = 0 (3.38)

Here θ stands for thrombin concentration, P0 denotes the concentration of platelets
at t = 0 and α > 0 is an activation rate. Arguing as for (3.4) under assumption
(3.3), equations (3.37), (3.38) yield at once:

pa(t) =

{
0 if t < TM

P0(1− e−αC) if t > TM
(3.39)

Upon activation, platelets are known to release a number of growth factors with
mitotic and chemotactic effects, including members of the VEGF, FGF, PADG,
and TGF families (cf. [1], [10], [13], [47]). We next describe the time and space
dynamics of one such factor g(x, t) irrespective of its specific biological function which
will be considered later on.

To that end we consider our fracture gap Ω = (−L
2 , L

2 ) with L > 0 and assume
that g(x, t) is released by activated platelets in Ω and freely diffuses within and
without Ω. More precisely, we postulate that g satisfies:

∂g

∂t
−D

∂2g

∂x2
+ βg = σ(t)pa(t)χΩ for −∞ < x < ∞, t > 0 (3.40)

g(x, 0) = 0 for −∞ < x < ∞ (3.41)

where χΩ(x) = 1 when x lies in Ω, and χΩ(x) = 0 otherwise. Here D and β denote
respectively the diffusion and half-life coefficients of g(x, t). Concerning σ(t), we
assume it to be of the form:

σ(t) =

{
σ > 0 for TM < t < T ∗

0 otherwise
(3.42)

for some constant σ and some T ∗ > 0. The reason for that choice is that platelets
can only release such substances as they had initially stored, since they lack nuclei
and hence the possibility to synthesize new molecules.

The linear problem (3.40), (3.41) can be explicitly solved by classical methods to
give:

g(x, t) =

∫ t

0

(∫ L/2

−L/2

G(x− y, t− s)dy

)
e−β(t−s)σ(s)pa(s)ds (3.43)

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 65

where

G(x, t) =
e−

x2

4Dt

√
4πDt

for −∞ < x < ∞, t > 0. (3.44)

Recalling (3.39), for times TM < t < T ∗ (3.43) yields:

g (x, t) = σP0

(
1− e−αC

) ∫ t

TM

(∫ L/2

−L/2

G (x− y, t− s) dy

)
e−β(t−s)ds (3.45)

whereas for t > T ∗:

g (x, t) = σP0

(
1− e−αC

) ∫ T∗

TM

(∫ L/2

−L/2

G (x− y, t− s) dy

)
e−β(t−s)ds (3.46)

3.3.4 Recruitment of mesenchymal cells (MSC). Formation of
early fibrin-collagen callus (EFC)

Having described the dynamics of a generic growth factor g (x, t) in our previous
Section, we now analyze the chemotactic effects induced by growth factors in the
population of MSCs and their impact on remodelling the fibrin net to produce a first
fibrin-collagen template. We assume that a baseline population of MSC is located at
the periosteum, near the fracture border and that when stimulated by g (x, t) such
population moves inwards into the fracture gap. When doing so, MSCs eventually
become precursors of various cell types instrumental to achieve bone repair. More
precisely, let us denote by m (x, t) the concentration of MSCs at a point x at a time
t. We postulate that m satisfies:

∂m

∂t
+

∂

∂x

(
µm

∂g

∂x

)
= 0 ; −∞ < x < ∞, t > TM (3.47)

m (x, TM ) = 0 for − L

2
≤ x ≤ L

2
, (3.48)

m (x, TM ) = m0 for x ≤ −L

2
, x ≥ L

2
. (3.49)

where µ > 0 is a chemotactic coefficient andm0 is the concentration of MSC available
at the periosteum reservoir. Notice that both diffusive and mitotic effects have been
ignored in (3.47). Using (3.45), (3.46) an explicit formula can be obtained for ∂g

∂x .
To this end, we note that on setting:

erf(u) =
2√
π

∫ u

0

e−x2

dx (3.50)

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 66

one has that

I (x, t) =

∫ L/2

−L/2

e−
(x−y)2

4Dt

√
4πDt

dy

=
1

2

(
erf

(√
Tc

t
(1/2 + x/L)

)
+ erf

(√
Tc

t
(1/2− x/L)

))
(3.51)

where

Tc =
L2

4D
. (3.52)

Then there holds:

∂I

∂x
(x, t) =

1√
π

√
Tc

t

1

L

(
e−

Tc
t (1/2+x/L)2 − e−

Tc
t (1/2−x/L)2

)
.

Setting P̂0 = σP0

(
1− e−αC

)
it follows that for TM < t < T ∗:

∂g

∂x
(x, t) = P̂0

∫ t

TM

∂I

∂x
(x, t− s) e−β(t−s)ds

=
P̂0√
πL

∫ t

TM

√
Tc

t− s

(
e−

Tc
4(t−s)

(1+2x/L)2 − e−
Tc

4(t−s)
(1−2x/L)2

)
e−β(t−s)ds . (3.53)

At this juncture it is convenient to make use of non-dimensional variables given by:

u =
2x

L
, τ =

t− s

Tc
, t̄ =

t− TM

Tc
(3.54)

Upon introducing (3.54) for times t̄ such that 0 < t̄ < T∗−TM

Tc
(3.53) can be recast

in the form:

∂g

∂x
(u, t̄) =

P̂0Tc√
πL

∫ t̄

0

(
e
−
(

1+u
2
√

τ

)2

− e
−
(

1−u
2
√

τ

)2
)

e−Tcβτ

√
τ

dτ (3.55)

When t̄ ≪ 1 so is τ , and the integrand in (3.55) is almost zero except for values close
to interval boundaries u = ±1. For instance, when t̄ = 0.01 a plot for ∂g

∂x (u, t̄) is
provided in figure 3.6 below.

When t̄ = O (1), that is for times t such that t−TM = O (Tc), the profile for
∂g
∂x (u, t̄)

looks as depicted in Figure 3.7. Such a profile is mantained until it begins to decrease
for t > T ∗. Therefore the main contributions to ∂g

∂x are achieved for TM < t < T ∗,

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 67

Figure 3.6: The form of ∂g
∂x for t̄ = 0.01. All parameters involved are set equal to

one for convenience.

Figure 3.7: A plot of ∂g
∂x for t̄ = 1. All parameters involved are set equal to one for

convenience.

Figure 3.8: A plot of ∂g
∂x (−1, t̄) as a function of t̄, showing quick stabilization to a

plateau-like profile. Parameter values as in Figures 3.6 and 3.7.

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 68

and are provided by values of u near to the boundary values u = ±1. For instance,
on writing:

∂g

∂x
(−1, t̄) =

P̂0Tc√
πL

∫ t̄

0

(
1− exp

(
−1

τ

))
e−Tcβτ

√
τ

dτ (3.56)

one may obtain a plot for ∂g
∂x (−1, t̄) as a function of t̄ that looks as in Figure 3.8:

Note that ∂g
∂x (−1, t̄) quickly saturates to a constant value. Actually, one has that:

∫ t̄

0

(
1− e−

1
τ

) e−Tcβτ

√
τ

dτ ≤ 1√
Tcβ

∫ ∞

0

(
1− e−

Tcβ
v

) e−v

√
v
dv

≤ 1√
Tcβ

∫ ∞

0

e−v

√
v
dv =

√
π

Tcβ

Hence:

µ
∂g

∂x
(−1, t̄) ≤ µP̂0

√
Tc/β

L
≡ µa , (3.57)

and a similar result can be derived at u = 1. We are now ready to estimate the time
it takes for a front of MSCs to reach the center of the gap at x = 0. Once front x (t)
has entered into a linear regime (which in particular occurs for times t− TM ≈ 3Tc;
see Figure 3.8), its equation reads:

dx

dt
= µ

∂g

∂x
(3.58)

x (TM ) = −L/2

Since ∂g
∂x ≈ −au = −2ax

L for 0 < t̄ < T∗−TM

Tc
and −L

2 < x < 0, (cf. 3.57) we actually
obtain that the differential equation in (3.58) can be approximated by:

dx

dt
= −2µa

L
x (3.59)

whence

x (t) = −L

2
e−(

2µa
L (t−TM )) (3.60)

On its turn, the concentration of MSC can be estimated from (3.59) and (3.47),
which together yield:

∂m

∂t
− 2aµx

L

∂m

∂x
=

2aµ

L
m (3.61)

Thus, if x (t) denotes the characteristic curve of (3.61) starting from x0 ≤ −L
2 , we

obtain from (3.61) that:

m (x (t) , t) = m0e
2µa
L (t−TM ) (3.62)

whereas m (x (t) , t) = 0 if x0 > −L
2 . We may now estimate the time TF of formation

of EFC as follows. Due to the approximations made to derive (3.60) the function

Universidad Complutense de Madrid



3.3. A MATHEMATICAL MODEL FOR THE FORMATION OF A FIBRIN-COLLAGEN CALLUS
TEMPLATE. 69

x(t) defined there will never reach x = 0 in a finite time. We will thus say that EFC
has been formed when x(t) reaches instead a value −ε with 0 < ε ≪ 1 arbitrarily
selected. We then may compute TF from the relation:

−L

2
e−

2µa
L (TF−TM ) = −ε

which gives at once:

TF = TM +

√
βDL

µP̂0

log

(
L

2ε

)
, (3.63)

provided that T ∗ > TF (see (3.55)).
We conclude this section by remarking about the negligible impact of diffusion

in the process under consideration. To avoid technicalities, we just show here that
adding diffusion in (3.47) results in a small correction in the inward motion of MSC
front. More precisely, let us replace (3.47) by

∂m

∂t
+

∂

∂x

(
µm

∂g

∂x
−Dm

∂m

∂x

)
= 0 ; −∞ < x < ∞ , t > TM (3.64)

where Dm > 0 with initial conditions (3.48) (3.49). Using again the approximation
∂g
∂x ≈ − 2a

L x in (3.64) the resulting equation can be written in the form

∂m

∂t
− 2µax

L

∂m

∂x
− 2µa

L
m−Dm

∂2m

∂x2
= 0 (3.65)

We first estimate the characteristic time τ trϵ needed to travel a distance L
2 under

transport effects. The time to fill an infinitesimal interval of the form [x− dx, x] is
given by:

dx

−2µax/L
= dt .

Integrating this expression between L/2 and a characteristic size epsilon (ϵ) (which
can be interpreted as an intercellular space) we find that τ trϵ satisfies the following
relationship: ∫ τtr

ϵ

0

dt =

∫ L
2

ϵ

L

2µax
dx

and then, we have that:

τ trϵ =
L

2µa
log

(
L

2ϵ

)
Analogously, we estimate the characteristic time τdiffϵ for the gap filling due to a
diffusion process:

τdiffϵ ∼
(
L
2 − ϵ

)2
Dm

Comparing both values of characteristic times we obtain:

τ tr

τdiff
=

LDm

2µa

log
(
L
2ϵ

)(
L
2 − ϵ

)2 ≈ LDm

2µa

log
(
L
2ϵ

)(
L
2

)2 =
2Dm

µaL
log

(
L

2ϵ

)

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 70

Assuming a value of ϵ of L
100 , then we can estimate the quotient (let us recall that

L = 2 mm):

τ tr

τdiff
=

2Dm

µaL
log

(
L

2ϵ

)
=

2 ∗ 10−6

0.5 ∗ 10−2 ∗ 1.25 ∗ 0.2
log(50) ≈ 6 · 10−3

which confirms the low impact of diffusion in the process under consideration, show-
ing that the transport effects occur in a much shorter characteristic time than diffu-
sive ones.

3.4 A sequence of times

We now briefly discuss on the comparative values of the characteristic times obtained
in our previous sections. To begin with, we observe that a correct timing of the overall
process requires that

TM < Tf < Tc < TF < T ∗ (3.66)

where TM , Tf , Tc, TF , and T ∗ have been respectively defined in (3.3), (3.21), (3.52),
(3.63) and (3.42). For simplicity, we take TM = 0 in the sequel, since hematoma
formation is known to begin shortly after fracture. On the other extreme, it is
reasonable to select T ∗ of the order of the average platelet life-span, so that we may
take:

T ∗ ≈ 10 days (3.67)

(cf. [33]). To estimate the values of the remaining times in (3.66) is in principle
feasible by means of the corresponding formulas obtained. However, while some of
the parameters appearing there can be reasonably guessed (at least from in vitro
experiments), little can be said about parameters as α and σ appearing in (3.37)
and (3.40) respectively. We next describe the type of information on the sequence
of times in (3.66) that can be obtained under current uncertainty constraints.

To begin with, consider the value obtained in (3.21) for Tf . In the earliest de-
tectable stages in fibrin clot formation, when fragmentation and deactivation effects
are ignored, the value of Tf can be compared with the prothrombin time Tp, a pa-
rameter routinely measured in blood tests [4]. Tp is measured in decalcified plasma
where platelets have been removed, and a typical value obtained is Tp ≈ 10 seconds.
If we compare the fastest estimate for Tf obtained in Section 3.2 (see (3.21)), we
should have:

Tf =
1

8FMkp (1− e−kgC)
≈ 10 sec (3.68)

FM in (3.68) is commonly estimated in blood tests. A typical value for FM is

FM ≈ 3 gr/l = 9000 nM (3.69)

(see [26]). Fibrin polymerization rates seem to have been reported only in vitro ([26]
and references therein). For instance, we may take

kp = 4× 10−3 (nM ·min)−1 (3.70)

Universidad Complutense de Madrid



3.4. A SEQUENCE OF TIMES 71

(cf. [58]). One then readily sees that 8FMkp ≈ 4.8 (sec)−1. Then estimate (3.68)
follows provided that

1− e−kgC =
1

48
(3.71)

which holds provided that kgC ≈ 0, 0211. As a matter of fact, values for kg have
been reported ([58], [26]). For instance, it has been suggested that

kg = 3× 10−4 (nM ·min)−1 (3.72)

(cf. [58]). From (3.71) and (3.72) we obtain that under such assumptions

C ≈ 104

3
0.0211 nM ·min = 0.47FM sec ≈ FM

2
sec . (3.73)

It has to be stressed that, as remarked in Section 2, the time to develop a mature
fibrin clot is considerably larger than the value provided in (3.68). Indeed, to keep
track of the actual unfolding of a fibrin net, quantities other than the prothrombin
time are measured in clinical practice [46]. Examples are the clotting time, CT,
defined as the time required for blood to form a initial clot from a given volume
of blood in a glass tube, which usually ranges between 5 to 15 minutes, and the
retraction time RT, defined as the time taken for a blood sample to become solid
enough so that it could separate from the walls of a crystal vessel where coagulation
is taking place, and that ranges between 1 to 2 hours. As a matter of fact, times as
those reported for CT or RT can be obtained from our formulas (3.21) or (3.25) by
selecting C and kb in a a suitable manner there.

Consider now Tc given in (3.52) through the relation

Tc =
L2

4D
(3.74)

where D is the diffusion coefficient of the growth factor considered. Since not much
seems to be known about in vivo values for D, we may now try on (3.74) a typical
diffusivity value, say D = 10−6 cm2/sec. Taking now L = 2 mm (see Appendix A)
we obtain:

Tc = 104 sec ≈ 3 hours (3.75)

Notice that by (3.68) and (3.75), Tf ≪ Tc. Moreover, a slower diffusivity will result
in an increased value of Tc according to (3.74).

Let us now remark on TF , the time required for early callus formation, for which

we have obtained estimate (3.63) in Section 3.4. For definiteness, let us take ε = e−N

2
for some integer N there. Or taking N such that N ≫ log(L), we see that

TF ≈
√
βDLN

µσP0 (1− e−αC)
(3.76)

Some of the parameters in (3.76) are known within reasonable ranges. For instance,
values for half-lives of growth factors have been reported in the literature, although

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 72

they considerably differ from one reference to another. For definiteness we take here
β = 6 · 10−4 (min)−1 [3]. Keeping to our previous choices L = 2 mm and D = 10−6

cm2/sec, and selecting N = 50√
6
we deduce that

√
βDLN ≈ 10−4 cm2/sec. Platelet

counts are commonly made in blood tests, with values ranging between 1.5 · 105 to
40 · 105 units per cubic millimeter, and mean platelet weight has been estimated to
be about 10 picograms [34]. We therefore may take P0 = 2gr/liter as a reference
value. On the other hand, the chemotactic coefficient µ has been measured for
neutrophils in some experimental settings [56] where values as 400 cm2/(M · sec)
have been reported. Just for comparison purposes, we may consider here growth
factor concentrations in the previous estimate for µ so that µP0 = 10−5 cm2/sec. In
this manner we will obtain: √

βDLN

µP0
≃ 10 (3.77)

Concerning the remaining parameters in (3.76) it is natural to assume σ = O
(

1
T∗

)
where T ∗ is given in (3.67). On taking σ = 1

T∗ and selecting then the conversion
factor 1− e−αC ≃ 1

50 , we would eventually arrive at

TF ≃ T ∗

5
≃ 2 days (3.78)

Needless to say, our previous estimates have to be considered as very preliminary.
For instance, for times of the order of CT or RT, additional characteristic times of
the various coagulation and fibrinolytic networks involved need to be investigated.
Actually, CT is used as part of hemophilic test, so that the dynamics of individual
coagulation factors as FVIII need to be accounted for there. On the other hand, once
the fibrin clot has been formed (which begins to occur quite fast; see [5]) the half life
of growth factors can be extended due to their anchorage in such a net ([49], [60]).
Actually, (3.76) suggests that an increase in that value from β to 4β would result in
doubling the corresponding value in (3.76). This effect, however, can be somehow
compensated by the slowdown in diffusivity D as the medium becomes denser.

3.5 Concluding remarks

In this work we have used mathematical models to explore the early stages of bone
fracture repair. Bone health is essential to ensure high quality of human life, and
its significance becomes even more important in aging societies, as many developed
countries are. Human bone has an impressive capability for self-repair. In fact,
human skeleton is constantly being remodelled, and as we have already noticed at
the Introduction, bone microfractures induced by exercise are routinely self-healed
without being ever noticed [6]. It has been observed that such repairing mechanisms
repeat many processes present at the developmental program that drives ossification
during the embryonic stage [18], [19]. Moreover, bone repair has some common
features with other tissue repairing processes as wound healing [32], [52]. It follows
from these remarks that each improvement in our understanding of the process of

Universidad Complutense de Madrid



3.5. CONCLUDING REMARKS 73

bone self-repair may have significant consequences, and it is widely expected that
quantitative models of bone healing may provide valid contributions to that end.

However, as in many other biological processes, when faced with the goal of pro-
viding mathematical models for bone repair, one has to confront a significant obsta-
cle. More precisely, while many of the key agents in the process under consideration
have been identified, little is known concerning the precise manner whereby these
agents coordinate themselves to achieve bone repair. For instance, only partial infor-
mation seems to be available on the detailed network hierarchy of the corresponding
biochemical cascades, and precise information on the order and rate constants of
the chemical reactions involved in cell signalling during that process remain largely
unknown. This leaves few options except selecting ad hoc such parameters when it
comes to carry out simulations to recover observed behaviours. However, it is well
known that a given dynamics can be recovered from different sets of parameter val-
ues, so that the insight gained through extensive parameter fitting has to be carefully
weighted ([29], [22]).

Bearing these remarks in mind, we have aimed in this work at discussing simple
models that could however provide some information into the specific issue addressed,
namely the timing of events taking place at the early stages of bone fracture repair.
Specifically, we restrain our attention to processes unfolding from bone break up
until the formation of a first fibrous callus that, upon suitable remodelling , will
lead to successful bone repair. While comparatively simple when compared with
later stages (for instance, mechanical effects do not play yet so relevant a role as in
posterior bone remodelling ) , we already deal here with some of the main processes
that will be at work later on, as chemotactic recruitment of MSCs (mesenchymal
stem cells). MSCs will thus divide and differentiate to produce cell lineages whose
interaction with previous elastic scaffolds will eventually lead to complete fracture
healing. Actually , the formation of the first (and crucial) such scaffold (a fibrin net)
is addressed in detail here , and its onset is kept track of and characterized as the
unfolding of a sol-gel phase transition.

Concerning the simplicity of our models, a few words are in order. When dealing
with any of the processes considered, we have attempted to focus on the main fea-
tures needed to derive the sought-for information, namely the time it takes for such
a process to develop. When doing so, we have considered effective processes, where
in particular various biochemical cascades have been lumped into a single equation,
rather than exploring the intricacies of each aspect involved. Our approach is illus-
trated by our discussion on thrombin-induced fibrin polymerization in Section 3.1,
where further remarks on the issues addressed can be found.

As a consequence of our work, some estimates for what we retain to be the main
characteristic times of the earliest stages considered have been obtained. These have
been summarized and estimated in Section 3 and 4 following a brief survey on the
main underlying biological aspects which makes the content of Section 2 . We have
discussed in section 4 how meaningful estimates for such characteristic times can be
provided, even in the presence of large parameter uncertainty, when our model is
compared with experimental data. However, additional work is required to improve

Universidad Complutense de Madrid



CHAPTER 3. EARLY STAGES OF BONE FRACTURE HEALING 74

the form of the estimates therein derived. Particularly, we have selected for that
purpose the induction of thin fractures in rats, for which a number of experiments
have been performed (for the corresponding experimental set up, see Appendix A).
Finally, while keeping models as simple as possible , we have succinctly discussed
on a number of possible extensions of the results obtained, as for instance the case
where thrombin sources are spatially distributed in Section 3.1, or a justification
of the negligible nature of diffusion aspects in the chemotactic recruitment of MSCs
during early callus formation in Section 3.4. Other aspects, however, will require of a
more detailed analysis. A particularly relevant example of the latter is the evaluation
of the robustness and redundancy properties induced by the simultaneous action of
several growth factors whith partially overlapping functions. We intend to address
this and other related issues in future work.

Universidad Complutense de Madrid



Bibliography

[1] Anitua, E., Andia, I., Ardanza, B., Nurden, P., Nurden, A.: Antologous
platelets as a source of proteins for healing and tissue regeneration. Thromb.
Haem. 91(1), 4–15 (2004)

[2] Ayati, B., Edwards, C., Webb, G., Wilswo, J.: A mathematical model of bone
remodelling dynamics for normal bone cell populations and myeloma bone dis-
ease. Biology Direct 5 (2010). http://www.biology-direct.com/content/5/1/28

[3] Beenken, A., Mohammadi, B.: The fgf family: Biology, pathophysiology and
therapy. Nat. Rev Drug. Discover. 8(3), 235–253 (2009)

[4] Van den Besselaaar, A.M.: Precision and accuracy of the international normal-
ized ratio in oral anticoagulant control. Haemostasis. 26(2s), 248–265 (1996)

[5] Bolander, M.: Regulation of fracture repair by growth factors. Proc Soc Exp
Biol Med 200, 165–170 (1992)

[6] Boyde, A.: The real response of bone to exercise. J. Anat 203(2), 173–189
(2003)

[7] Brass, L.: Thrombin and platelet activation. Chest 124(3), 185–255 (2003)

[8] Butenas, S., Vant Veer, C., Mann, K.: Normal thrombin generation. Blood
94(7), 2169–2178 (1999)

[9] Chung, R., Cool, J., Scherer, M., Foster, B., Xian, C.: Roles of neutrophil-
mediated inflammatory response in the bony repair of injured growth plate
cartilage in young rats. J. Leukoc. Biol 80, 1272–1280 (2006)

[10] Cicha, I., Garlichs, C., Daniel, W., Goppelt-Struebe, M.: Activated human
platelets release connective tissue growth factor. Throm. Haem. 91(4), 755–760
(2004)

75



BIBLIOGRAPHY 76

[11] Colnot, C., Huang, S., Helms, J.: Analyzing the cellular contribution of
bone marrow to fracture healing using bone marrow transplantation in mice.
Biochem. Biophys Res Commun 350, 557–561 (2006)

[12] Danforth, C., Orfeo, T., Mann, K., Brummel?Ziedins, K., Everse, S.: The
impact of uncertainty in a blood coagulation model. Math. Med. Biol. 4, 323–
336 (2009)

[13] Dimitriou, R., Tsiridis, E., Giannoudis, P.: Current concepts of molecular as-
pects of bone healing. Injury Int. J. Care. 36, 1392–1404 (2005)

[14] Eghbali-Fatourechi, J., Lamsam, D., Fraser, D., Nagel, B., Riggs, S., Khosla:
Circulating osteoblast-lineage cells in humans. N Engl J Med. 352, 1959–1966
(2005)

[15] Einhorn, T.: The cell and molecular biology of fracture healing.
Clin.Orthop.Relat. Res 355(S), 7–21 (1998)

[16] Fasano, A., Herrero, M., Lopez, J., Medina, E.: On the dynamics of the growth
plate in primary ossification. Journal Theor. Biol 265, 543–553 (2010). doi:
10.1016/j.jtbi.2010.05.030

[17] Fasano, A., Santos, R., Sequeira, A.: Blood coagulation: a puzzle for biolo-
gists, a maze for mathematicians. In: Modelling Physiological Flows, pp. 44–76.
Springer Verlag (2011)

[18] Ferguson, C., Miclau, T., Alpern, E., Helms, J.: Does adult fracture repair
recapitulate embryonic skeletal formation? Mech. Dev. 87, 57–66 (1999)

[19] Ferguson, C., Miclau, T., Hu, D., Alpern, E., Helms, J.: Common molecular
pathways in skeletal morphogenesis and repair. Ann. New York Acad. Sci. 857,
33–42 (1988)

[20] Gerstenfeld, L., Alkhiary, Y., Krall, E., Nicholls, F., Stapleton, S., Fitch, J.,
Bauer, M., Kayal, R., Graves, D., Jepsen, K., Einhorn, T.: Three-dimensional
reconstruction of fracture callus morphogenesis. J HistochemCytochem 54,
1215–1228 (2006)

[21] Gerstenfeld, L., Cullinane, D., Barnes G.L.and Graves, D., Einhorn, T.: Frac-
ture healing as a post-natal developmental process: molecular, spatial, and
temporal aspects of its regulation. J Cell Biochem 88, 873–884 (2003)

[22] Ginzburg, L.R., Jensen, C.X.J.: Rules of thumb for judging ecological theories.
Trends in ecology and evolution. 19(3), 121–126 (2004)

[23] Graham, J., Ayati, B., Holstein, S., Martin, J.: The role ofosteocytes in targeted
bone remodelling: a mathematical model. PLoSOne 22(8(5e63884)) (2013). doi:
10.1371/journal.pone.0063884

Universidad Complutense de Madrid



BIBLIOGRAPHY 77

[24] Groothuis, A., Duda, G., Wilson, C., Thompson, M., Hunter, M., Simon, P.,
Bail, H., van Scherpenzeel, K., Kasper, G.: Mechanical stimulation of the
pro-angiogenic capacity of human fracture haematoma: involvement of vegf
mechano-regulation. Bone 47, 438–444 (2010)

[25] Guria, G., Herrero, M., Zolobina, K.: A mathematical model of blood coagula-
tion induced by activation sources. Discrete and Continuous Dynamical Systems
25(1), 175–194 (2009)

[26] Guria, G.T., Herrero, M.A., Zlobina, K.E.: Ultrasound detection of externally
induced microthrombi cloud formation: a theoretical study. J. Eng. Math. 66,
293–310 (2010)

[27] Guria, K., Gagarina, A., Guria, G.: Instabilities in fibrinolytic regulatory sys-
tem. theoretical analysis of blow-up phenomena. J. Theor. Biology 304, 27–38
(2012)

[28] Guy, R., Fogelson, A., Keener, J.: Fibrin gel formation in a shear flow. Math.
Med. Biol. 24, 111–130 (2007)

[29] Hemker, H.C., Kerdelo, S., Kremers, R.M.W.: Is there value in kinetic modelling
of thrombin generation? J. Thromb. Haemost. 10, 1470–1477 (2012)

[30] Hollinger, J., Onikepe, A., MacKrell, J., Einhorn, T., Bradica, G., Lynch, S.,
Hart, C.: Accelerated fracture healing in the geriatric, osteoporotic rat with
recombinant human platelet-derived growth factor-bb and an injectable beta-
tricalcium phosphate/collagen matrix. J Orthop Res 26, 83–90 (2008)

[31] Hughes, S., Hicks, D., OKeefe, R., Hurwitz, S., Crabb, I., Krasinskas, A.,
Loveys, L., Puzas, J., Rosier, R.: Shared phenotypic expression of osteoblasts
and chondrocytes in fracture callus. J Bone Miner Res 10, 533–544 (1995)

[32] Javierre, E., Vermolen, F., Vuik, C., Van der Zwaag, S.: A mathematical anal-
ysis of physiological and morphological aspects of wound closure. J. Math. Biol
59, 605–633 (2009). doi: 10.1007/s 00285-008-0242-7

[33] Kameth, S., Blann, A.D., Lip, G.Y.H.: Platelet activation: assesment and quan-
tificaion. European Hearth J. 22, 1561–1571 (2001)

[34] Kiem, J., Borberg, H., Iyengar, G.V., et al.: Water content of normal human
platelets and measurements of their concentrations og cu, fe, k and zn by neutron
activation analysis. Clin. Chem. 21(5), 705–710 (1979)

[35] Komarova, S., Smith, R., Dixon, J., Sims, S., Wahl, L.: Mathematical model
predicts a critical role for osteoclast autocrine regulation in the control of bone
remodelling. Bone 33, 206–215 (2003). doi: 10.1016/S8756-3282(03)00157-1

[36] Malizos, N., Papatheodorou, L.: The healing potential of the periosteum. molec-
ular aspects. Injury 36(S3), S13–S19 (2005)

Universidad Complutense de Madrid



BIBLIOGRAPHY 78

[37] Mann, K., Butenas, S., Brummel, K.: The dynamics of thrombin formation.
Arteriosclerosis Thrombosis and Vascular Biology 23, 17–25 (2003)

[38] Mc Dougall, S., Dallon, J., Sherratt, J., Maini, P.: Fibroblast migration and
collagen deposition during termal wound healing: mathematical modelling
and clinical implications. Phil. Trans. Roy. Soc 364, 1385–1405 (2006). doi:
10.1098/rsta.2006.1773

[39] McKibbin, B.: The biology of fracture healing in long bones. J. Bone Joint
Surg.B 2(60), 151–162 (1978)

[40] Moore, S., Saidel, G., Knothe, U., Knothe Tate, M.: Mechanistic , mathematical
model to predict the dynamics of bone defects via mechanical feedback and
mediation of biochemical factors. PLoS Comp. Biol 10(6e1003604) (2014). doi:
10.1371/journal.pcbi.1003604

[41] Moreo, P., Garcia Aznar, J., Doblare, M.: Bone ingrowth on the surface of
endosseous implants. part i: mathematical model. J. Theor. Biol 260, 1–12
(2009)

[42] Murphy, K., Hall, C., Maini, P., Mc Cue, S., Sean Mc Elwain, D.: A fibrocon-
tractive mechanochemical model of wound closure incorporating realistic growth
factor kinetics. Bull. Math. Biol 74, 1143–1170 (2012). doi: 10.1007/s11538-011-
9712-y

[43] Orfeo, T., Butenas, S., Brummel-Ziedins, K., Mann, K.: The tissue factor
requirement in blood coagulation. J. Biol. Chemistry 280(52), 42,887–42,896
(2005)

[44] Pivonka, P., Dunstan, C.: Role of mathematical modeling in bone fracture
healing. BioKEyReports 1(221) (2012). doi: 10.1038/bonekey 2012.221

[45] Prokharau, P., Vermolen, F., Garcia-Aznar, J.: Model for direct bone apposition
on a pre-existing surface during peri implant osteointegration. J. Theor. Biol
304, 131–142 (2012)

[46] Rodak, B.F., Fritsma, G.A., Keohane, E.: Hematology: Clinical Principles and
Applications., 4 edn. Saunders, USA. (2011)

[47] Rozman, P., Boita, Z.: Use of platelet growth factor in treating wounds and
soft-tissue injuries. Acta. Dermat. 16(4), 156–165 (2007)

[48] Rumi, M., Deol, G., Singapuri, K., V.D., P.J.: The origin of osteo- progenitor
cells responsible for heterotopic ossification following hip surgery: an animal
model in the rabbit. J Orthop Res. 23, 34–40 (2005)

[49] Sahni, A., Francis, C.W.: Vascular endothelial growth factor binds to fibrinogen
and fibrin and stimulates endothelial cell proliferation. Blood 96, 3772–3778
(2000)

Universidad Complutense de Madrid



BIBLIOGRAPHY 79

[50] Schmidt-Bleek, K., Schell, H., Schulz, N., Hoff, P., Perka, C., Buttgereit, F.,
Volk, H., Lienau, J., Duda, G.: Inflammatory phase of bone healing initiates
the regenerative healing cascade. Cell Tissue Res 347, 567–573 (2012)

[51] Shapiro, F.: Bone development and its relation to fracture repair. the role of
mesenchymal osteoblasts and surface osteoblasts. Eur Cells Mater. 15, 53–76
(2008)

[52] Sonnemann, K.J., Bement, W.M.: Wound repair: toward understanding and
integration of single cell and multicellular wound responses. Ann. Rep. Cell and
Develop Biol. 27, 237–263 (2011)

[53] Stockmayer, W.: Theory of molecular size distribution and gel formation in
branched-chain polymers. J.Chem. Physics 11(2), 45–55 (1943)

[54] Street, J., Winter, D., Wang, J., Wakai, A., McGuinness, A., Redmond, H.: Is
human fracture hematoma inherently angiogenic? Clin Orthop 378, 224–237
(2000)

[55] Takahara, M., Naruse, T., Takagi, M., Orui, H., Ogino, T.: Matrix
metalloproteinase-9 expression, tartrate-resistant acid phosphatase activity, and
dna fragmentation in vascular and cellular invasion into cartilage preceding pri-
mary endochondral ossification in long bones. J Orthop Res. 22, 1050–1057
(2004)

[56] Tranquillo, R.T., Zigmund, S.H., Lauffenburger, D.A.: Measurement of the
chemotaxis coefficient for human neutrophils in the under-agarose migration
assay. Cell Motility and Cytoskeleton 11, 1–15 (1988)

[57] Wagenwood, R., Hemker, P.W., Hemker, H.C.: The limits of simulation of the
clotting system. J. Thromb. Haemost. 4, 1331–1338 (2006)

[58] Weisel, J.W., Veklich, Y., Gorkun, O.: The sequence of cleavage of fibropeptides
from fibrinogen is important for protofibrin formation and enhacement of lateral
aggregation. J. Mol. Biol. 232, 285–297 (1993)

[59] Wiltzius, P., Dietler, G., Kanzig, W., Haberli, A., Straub, P.: Fibrin polymeriza-
tion studied by static and dynamic light scattering as a function of fibrinopeptide
a release. Biopolymers 21, 2205–2223 (1982)

[60] Yon, Y.R., Won, J.E., Jeon, E., et al.: Fibroblast growth factors: biology,
function and applications for tissue engineering. J. Tissue Eng. 1 (2010)

[61] Yoo, J., Johnstone, B.: The role of osteochondral progenitor cells in fracture
repair. Clin. Orthop. Relat. Res. 3(55), S73–S81 (1998)

Universidad Complutense de Madrid



BIBLIOGRAPHY 80

Universidad Complutense de Madrid



Appendix A
Details of computer simulations of
the model of BMU operation during
bone remodeling

A.1 Time evolution of the BMU

• Initialization: Simulations take place in a two-dimensional cellular automata
(CA) of 150 x 250 boxes. 125 rows of this matrix are initially occupied by bone
matrix and the rest by bone marrow. A row of OBL lies in the interface between
them. OCYs are initially located at fixed, regularly spaced positions within
the bone matrix. A random perturbation is then applied to their positions
to avoid any potential effect of symmetries on the simulations of the model.
At this point, OCYs and OBLs release signals that diffuse through the bone
matrix until an equilibrium is attained. This equilibrium is determined by the
balance between signal production and decay. In this moment, a region of OCY
apoptosis is randomly defined within a priori defined boundaries (Depth :=
maximum depth and Width := maximum width).

• Time evolution:

Simulations of the model occur in discrete time steps of length ∆t = 0.01. The
input for an iteration of the model corresponding to time t = t0 includes a set
of vectors (one for each of the cells coexisting in the model at that time) with
the coordinates of the box occupied by the cell, and its state variables (amount
of inhibitory molecules).

1.- OBL

The state variable considered in the model for the i-th OBL is the amount of
activation inhibitor dBi . The corresponding vector is:

81



APPENDIX A. DETAILS OF COMPUTER SIMULATIONS OF THE MODEL OF BMU
OPERATION DURING BONE REMODELING 82

vBi(t0) = {(x, y), dBi(t0)}

2.- OBA

For each OBA there are three different inhibitory molecules, namely cA, dA
and aA, blocking cell division, differentiation and apoptosis. The state vector
for the i-th OBA is:

vAi(t0) = {(x, y), {cAi(t0), dAi(t0), aAi(t0}}

3.- OPL

Differentiation of OPL into OCL is controlled by the inhibitor dP , so that the
state vector of the i-th OCP is given by:

vPi(t0) = {(x, y), dPi(t0)}

4.- OCL

Active osteoclast can die by apoptosis and dig into old bone. The state vector
for the i-th OCL includes the amount of apoptosis inhibitor aCi

and the amount
of old bone bCi at the box occupied by the OCL:

vCi
= {(x, y), {aCi

(t0), bCi
(t0)}}

Input also includes the values of signals S, T and R at each box of the CA.

The first step of each iteration consists in the definition of a system of ordinary
differential equations that describe the dynamics of the state variables (see
Section 2). It consists of NB(t0)+3×NA(t0)+NP (t0)+2×NC(t0) equations,
where NB, NA, NP and NC are the number of OBL, OBA, OCP and OCL at
time t0 respectively. The system is then solved numerically using the explicit
Runge-Kutta method implemented in the function NDSolve of Mathematica
(Wolfram Research, Inc., Mathematica, Version 9.0, Champaign, IL).

If neither the amount of bone at a box occupied by a OCL or some of the
inhibitors goes to zero during a time step the state of signals S, R and T
is actualized and the model moves to the next time step. Alternatively, the
model displays different behaviors depending on what variable (or variables)
have gone to zero:

1.- OBLs activation inhibitor (dB)

This corresponds to the activation of one or more OBLs. In this case, the acti-
vated OBL are removed from the list of cells of the model, and their positions
in the CA are occupied by newly formed OBAs.

2.- OBAs division inhibitor (cA)

Universidad Complutense de Madrid



A.1. TIME EVOLUTION OF THE BMU 83

OBAs form a continuous layer of cells that progresses into the cutting cone
created by OCL. By construction, only OBAs that are not located in boxed
adjacent to bone can divide. Moreover, OBA division can only occur if some of
the OBA are not attached to bone yet. In case of OBA division, the position
of OBAs is recalculated as summarized in Figure A1.

Figure A.1: Movement of OBAs after cell division. A) The layer of OBAs (in
orange) is adjacent to old bone (in grey). There are 7 available positions for new
OBAs marked with a thick black line. B) OBAs attached to the bone surface (in
red) can neither divide nor move. C) In case of cell division, one new OBA appears.
In the case depicted in the figure there are 6 OBAs after cell division that can move,
and 7 free positions, so the layer of OBAs cannot progress into the resorpted bone.
The new OBA is located in a cell adjacent to that of the dividing OBA. D) After a
new cell division there are 7 OBA free to move, so the front of OBA can occupy the
next set of available positions as shown in E. E) OBAs that reach a position adjacent
to bone (OBAs labelled with 1 and 7 in red) remain fixed and do not divide. F)
When all OBAs are attached to the bone surface cell division is no longer possible.
From this moment, OBAs start to secrete osteoid matrix, can undergo apoptosis and
move downwards.

3.- OBAs apoptosis and differentiation inhibitors (aA and dA)

Universidad Complutense de Madrid



APPENDIX A. DETAILS OF COMPUTER SIMULATIONS OF THE MODEL OF BMU
OPERATION DURING BONE REMODELING 84

During the bone formation stage, OBAs form a continuous layer of cells that
fills the cutting cone created by OCL with new osteoid matrix. Apoptosis of
OBAs and differentiation into OCY entail the movement of the OBA layer into
the empty cutting cone as described in Figure A.2. In case of differentiation,
the OBA is replaced by a new OCY in the same box.

Figure A.2: Movement of OBAs after apoptosis and differentiation. A) A
layer of 15 OBAs (in red) is adjacent to the bone matrix. There are 12 available
positions for the new layer of OBAs, marked with a thick black line. New positions
are bounded by the layers of OBAs and OBLs. B) In this moment, the number
of apoptotic (8 and 12, in orange) and differentiated (4, in blue) OBAs equals the
difference between the current number of OBAs (15) and the number of available
positions (12), so that the layer of OBAs can move. C) Boxes that were occupied
by OBAs in the previous stage are now filled with osteoid (yellow). A new OCY is
placed in the box where an OBA has differentiated.

4.- OCPs differentiation inhibitor (dP )

OCPs are located adjacent to the layer of OBLs. In case of activation, the
corresponding OCP is replaced by a newly formed OCL that moves to the
closed box within the bone matrix that is not occupied by any other cell.

5.- OCLs apoptosis inhibitor (aC)

OCLs whose apoptosis inhibitor disappears are removed from the set of cells
in the model for the next iteration.

6.- Bone removal by OCPs (bP )

OCLs remove old bone from the box where they are located according to
equation []. When the amount of bone goes to zero, the OCL moves to one of
the unoccupied adjacent boxes in the bone matrix. The OCL choses the box
with the highest amount of S signal that is farther from the OBL layer.

After cell fate choices, the set of cells and the amounts of signals S, T and R
are updated and the model starts a new iteration.

• End of the simulation:

Simulations end when the bone remodeling process is completed (old bone has
been replaced by new osteoid matrix) and a new equilibrium (for signals S and

Universidad Complutense de Madrid



A.1. TIME EVOLUTION OF THE BMU 85

T ) has been reached.

A.1.1 Pseudo-code guidelines

Input parameters

tmax := maximum duration of simulations

Depth := maximum depth of the apoptosis region

Width := maximum width of the apoptosis area

bone(x, y) := 1 if the box is occupied by bone at initial time and 0 other-
wise

position(OCY, 0) := initial position of OCYs

position(OBL, 0) := initial position of OBLs

position(OCP, 0) := initial position of OCPs

Structural parameters

αh
k := effect of external signals of type h on the dynamics of inhibitors in

cells of type k

DB := maximum amount of inhibitor dB that can accumulate inside a
OBL

dph(d) := diffusion profile of signal h as a function of distance d

γh := rate of decay of signal of type h

Variable parameters

{ak0, dk0, ck0} := amount of inhibitors in newly formed cells of type k

δA Distance of OCY inhibition of OBA differentiation

Qh
k := rate of production of signal of type h by cells of type k

Initialization

1 Input the bone matrix bone(x, y)

3 Input and the boxes occupied by OCY, OBL and OCL.

4 Diffusion and decay of signals S and T until equilibrium

5 Define a random region of OCY apoptosis

6 Set t := 0

7 REPEAT

8 Input the set of positions and state variables of all cells

9 Input the amount of signals S, T and R at each box of the CA

Universidad Complutense de Madrid



APPENDIX A. DETAILS OF COMPUTER SIMULATIONS OF THE MODEL OF BMU
OPERATION DURING BONE REMODELING 86

10 Write the system of ODEs that describe the dynamics of the state
variables as functions of the external signals

11 REPEAT

12 Numerical integration of the systems of ODEs for a time step ∆t

13 Diffusion and decay of signals S, T and R during the time step
∆t

14 set t = t+∆t

15 UNTIL some state variable goes to zero

16 If (dBi(x, y; t) ≤ 0) then remove the i-th OBL and add a new OBA
at box (x, y)

17 If (cAi(x, y; t) ≤ 0) then

18 add a new OBA at a box adjacent to (x, y)

19 recalculate OBAs positions

20 If (aAi(x, y; t) ≤ 0) then

21 remove the i-th OBA

22 recalculate OBAs positions

23 If (change in OBAs position) then fill old positions with osteoid

24 If (dAi(x, y; t) ≤ 0) then

25 remove the i-th OBA

26 add a new OCY at (x, y)

27 recalculate OBAs positions

28 If (change in OBAs position) then fill old positions with osteoid

30 If dPi(x, y; t) ≤ 0 then

31 remove the i-th OCP

32 add a new OCL at the box of the bone matrix closest to (x, y)

33 If aCi(x, y; t) ≤ 0 then remove the i-th OCL

34 If bBi(x, y; t) ≤ 0 then

35 set bone(x, y) = 0

36 move the OCL to a new position in the bone matrix

# We remark that several cell decisions can take place simultaneously

37 UNTIL the end of bone remodeling process

Universidad Complutense de Madrid



A.1. TIME EVOLUTION OF THE BMU 87

Parameter Description Value

(Depth,Width)
Maximum depth and widht of the

OCY apoptosis region
(80, 80)

dB0
Initial amount of differentiation
inhibitor in newly formed OBLs

3

DB
Maximum amount of differentiation

inhibitor in OBLs
4

(αS1

B , αS1

B )
OBL

structural parameters
(1, 1)

(cA0, aA0, dA0)
Initial amount of inhibitors in

newly formed OBA
(1.5, 5, 2)

(αT
A, α

S1

A , αS2

A , αS3

A , αS4

A )
OBA

structural parameters
(0.01, 0.6, 15, 1)

δA
Radio of OCY ihibition
of OBA differentiation

1

dP0
Initial amount of differentiation

inhibitor in OCP
0.5

(αR
P , α

S
P )

OCP
structural parameters

(3, 5)

aC0
Initial amount of apoptosis

inhibitor in OCL
1

(αS1

C , αS2

C , αS3

C , αS4

C )
OCL

structural parameters
(2, 2.5, 2, 6)

(dpS(d), dpT (d), dpR(d)
Diffusion profile of
signals S, T and R

dpS(d) = dpT (d) =
= dpR(d) = Max(0, 4− d)

(γS , γT , γR)
Decay rates of

signals S, T and R
(0.01, 0.01, 0.01)

(QS
Y , Q

T
Y )

Rates of production of
signals S and T by OCYs

(10, 1)

QT
B

Rate of production of
signal T by OBL

2

(QT
A, Q

R
A)

Rates of production of
signals T and R by OBAs

(1, 15)

Table A.1: Parameters

Universidad Complutense de Madrid



APPENDIX A. DETAILS OF COMPUTER SIMULATIONS OF THE MODEL OF BMU
OPERATION DURING BONE REMODELING 88

initialize 

structural 
parameters 

variable 
parameters 

initial CA 
configuration 

OCY 
apoptosis 

calculate 
signals 

equilibrium 
write the 
system  

of ODEs 

write the

numerical 
integration 

t = t + t 

all state  
variables 

positive?  

YES update 
signals 

YESS 

t =t

a

cells 
decisions 

NO NO 

OCP 
activation 

OBL 
activation 

OBA 
division 

OBA 
apoptosis 

OCL 
bone 

removal 
OBA 
diff 

ecisions de

OOCO PP
ctivatio

OBL
activationac

OBA
divisiond

OBA
apoptosis a

on

OOCO LL
bone 

removal 
OBA
diffffdifff

removal p p

new 
OCL OCL

remove 
OCL 

new 
OCY 

new 
OBA 

remove 
OBA 

update 
OCL position 

update 
bone 

update 
OBA positions ns 

OBA

u
b

re
OOCY

te
e

ac

OCL 
apoptosis 

update
OC tiCL positL

Figure A.3: Flowchart of the proposed model of BMU software as described
in the Appendix.

Universidad Complutense de Madrid



Appendix B
Experimental procedure

Sixteen 10-week-old male Sprague-Dawley rats were obtained from the animal fa-
cility building of the University of Oviedo. Bone fractures with 2.5 mm average
width were produced by the manual breakage using plate-bending devices, placed
at the distal 3 third of the right tibia. The animal study protocol was approved
by the local committee and conformed to Spanish animal protection laws. Rats
were anesthetized and sacrificed by cervical dislocation 3, 6, 12 and 24 hours af-
ter fracture (n=5). Tibiae were isolated and cut through the sagittal plane of the
metaphyseal fracture region into two equal sized parts and fixed by immersion in 4%
paraformaldehyde at 4 Co for 12 h, rinsed in PBS, decalcified in 10% EDTA (pH 7.0)
for 48 hours at 4 Co, dehydrated through a graded ethanol series, cleared in xylene
and embedded in paraffin. Sections were serially cut at a thickness of 5 m, mounted
on Superfrost Plus slides (Menzel-Glaser), stained with a trichromic method using
Weigert?s hematoxylin/alcian blue/picrofuchsin and viewed and photographed on a
Nikon Eclipse E400 microscope.

89


	Tesis Luis Echeverri Delgado
	Portada
	Contents
	Resumen
	Abstract
	Chapter 1.- Preliminaries
	Bibliography

	Chapter 2.- Bone remodelling
	Bibliography

	Chapter 3.- Early stages of Bone fracture healing
	Bibliography

	Appendix A.- Details of computer simulations ofthe model of BMU operation during bone remodeling
	Appendix B.- Experimental procedure