A bone remodeling model governed by cellular micromechanics and physiologically based pharmacokinetics

Journal Pre-proof A bone remodeling model governed by cellular micromechanics and physiologically based pharmacokinetics M.T. Bahia, M.B. Hecke, E.G.F. Mercuri, M.M. Pinheiro PII:

S1751-6161(19)31620-0

DOI:

https://doi.org/10.1016/j.jmbbm.2020.103657

Reference:

JMBBM 103657

To appear in:

Journal of the Mechanical Behavior of Biomedical Materials

Received Date: 23 October 2019 Revised Date:

11 January 2020

Accepted Date: 23 January 2020

Please cite this article as: Bahia, M.T., Hecke, M.B., Mercuri, E.G.F., Pinheiro, M.M., A bone remodeling model governed by cellular micromechanics and physiologically based pharmacokinetics, Journal of the Mechanical Behavior of Biomedical Materials (2020), doi: https://doi.org/10.1016/j.jmbbm.2020.103657. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

A bone remodeling model governed by cellular micromechanics and physiologically based pharmacokinetics M. T. Bahiaa,b,∗, M. B. Heckeb , E. G. F. Mercurib,c , M. M. Pinheirod a

Departamento de Mecˆanica, Instituto Federal de Santa Catarina, Rua Pav˜ao, 1377, Costa e Silva, CEP 89220-618, Joinville, Brasil b Grupo de Bioengenharia, Programa de P´os-Graduac¸a˜ o em M´etodos Num´ericos em Engenharia, Universidade Federal do Paran´a, Curitiba, Brasil c Departamento de Engenharia Ambiental, Universidade Federal do Paran´a, Curitiba, Brasil d Departamento de Medicina, Universidade Federal de S˜ao Paulo, S˜ao Paulo, Brasil

Abstract This study describes a mathematical model for bone remodeling that integrates the bone cells activities with the pharmacological dynamics for bone-seeking agents. The evolution of bone cells population involves the osteoblast-osteoclast signaling mediated by biochemical factors and receives both mechanical stimulus evaluated at the microscale and pharmacological regulation. A physiologically based pharmacokinetic model (PBPK) for bone-seeking agents was developed to provide the drug concentration on bone sites and feed the remodeling algorithm. The drug effect on bone was reproduced coupling three different strategies: modification of the RANKL expression, increase the osteoclast apoptosis and change in the rate of differentiation of preosteoblasts. Computational simulations were performed in the PBPK model considering different dosing regimens. A 3D finite element model of a proximal femur was generated and the simulation of the ∗

Corresponding author Email address: [email protected] (M. T. Bahia)

Preprint submitted to Journal of the Mechanical Behavior of Biomedical MaterialsJanuary 30, 2020

bone remodeling algorithm were implemented in Matlab. The results indicate that the proposed integrated model is able to capture adequately the expected adaptive behavior of bone subjected to mechanical and pharmacological stimulus. The model demonstrated to have potential for use as a platform to investigate therapies and may help in the study of new drugs for bone diseases. Keywords: bone cell dynamics, physiologically based pharmacokinetics (PBPK), continuum michomechanics, osteoporosis, mathematical modeling

1

1. Introduction

2

Bones are a living tissue capable of adapting mass and microstructure in re-

3

sponse to mechanical and biochemical stimuli in a physiological process called

4

bone remodeling. At the cellular level, remodeling involves primarily osteoblast-

5

mediated bone formation and osteoclast-mediated bone resorption. Osteocytes

6

acts as mechanotransducers receiving mechanical stimulus and producing bio-

7

chemical signaling agents (Fig. 1), such as nitric oxide (ON) and prostaglandin

8

E2 (PGE2 ), that will drive osteoblasts and osteoclast response (Rieger, 2011; Ross

9

et al., 2017; Qin and Wang, 2012). Bone diseases such as osteoporosis, osteope-

10

nia and Paget’s disease (Barkaoui, 2012) cause imbalances in bone remodeling.

11

In the case of osteoporosis, there is a degradation of the bone architecture and

12

structural properties that results in decreased bone quality and greater risk of

13

fractures (NIH, 2001).

14

The mortality rate among patients one year after osteoporotic fractures is

15

15 to 30 % due to complications such as infections, venous thrombosis, pressure

16

ulcers, or comorbid conditions, particularly cardiovascular diseases (Pinheiro,

17

2008; Pinheiro and Eis, 2010).According to Fortes et al. (2008) only 30 % of these 2

18

patients return to their daily activities and 11.6 % become totally dependent. In

19

Brazil, it is estimated that the cost of hospitalization due to osteoporotic fractures

20

in a private hospital is approximately US$ 12, 000 per patient (over 50 years of

21

age) (Ara´ujo et al., 2005). Osteoporosis represents a public health concern given

22

its high economic and social impact. With the increase in life expectancy, its

23

burden tends to increase.

Figure 1: Cartoon representing the bone remodeling process with the osteoblast-mediated bone formation and the osteoclast-mediated bone resorption driven by osteocytes. Osteocytes act in remodeling similarly to a master builder or a regent in an orchestra, signaling for osteoblasts and osteoclasts according to the mechanical/biochemical stimuli received. The cartoon shows a sleeping osteocyte, which is intended as a humorous content, as it is highly unlikely that an osteocyte would be slacking off amidst all the multiple signals it receives, be it direct strain, hydrostatic pressure, or interstitial fluid flow (Pastrama et al., 2018; Scheiner et al., 2016).Art: Miguel T. Bahia.

24

Given its financial and social relevance, the understanding of the processes

25

that lead to bone diseases and possibly to the development of new efficient ther3

26

apies are of paramount importance. Mathematical models of bone remodeling

27

have the ability to create virtual scenarios that allow the study of the dynamics

28

of bone cells populations, drug effects and different therapeutic strategies.

29

The approaches adopted to simulate bone remodeling numerically can be

30

quite different. Among them are mechanistics, phenomenologicals, multiscale,

31

and optimization methods, and hybrid techniques (Rieger, 2011; Doblar´e and

32

Garcıa, 2002; Garc´ıa-Aznar et al., 2005; Ganghoffer, 2012; Goda et al., 2016; Louna

33

et al., 2019). Mathematical models that consider the evolution of bone cells pop-

34

ulations or cell activities, include the pioneering works of Kroll (2000), which

35

analyzes the dual action of the parathyroid hormone (PTH) on bone formation

36

and resorption; Komarova et al. (2003), which models the autocrine and paracrine

37

interactions of osteoblasts and osteoclasts; and the model of Lemaire et al. (2004)

38

that explicitly expresses the RANK-RANKL-OPG pathway. Based on Lemaire

39

et al. (2004) model, the work of Pivonka et al. (2008) rewrites the signaling ex-

40

pressions with activation and repression functions from chemical kinetics (Hill

41

functions) and incorporates other signaling into the model. Another important

42

work is Scheiner et al. (2013) that combines the Pivonka et al. (2008) model with

43

the continuum micromechanics and incorporates the mechanical stimulus evalu-

44

ated in the microscale and use the Mori-Tanaka homogenization technique (Mori

45

and Tanaka, 1973).

46

Numerical simulations with pharmacological models represent a whole new

47

avenue of possibilities for bone remodeling studies. Parameters such as dose, fre-

48

quency, route of administration and effectiveness of a medicine can be included

49

in the model, making the simulation quite relevant in terms of treatment, preven-

50

tion and cure of bone diseases. The application of the pharmacological modeling 4

51

involves (i) the development of a pharmacokinetic model that simulates the tem-

52

poral evolution of the drug concentration in the body or a specific organ, (ii) the

53

development of a pharmacodynamic model that simulates the drug’s action in

54

the body or organ. In other words, how the drugs interacts with the body, and

55

how the body reacts to the drug. In both cases, the models should be able to

56

reasonably simplify a complex physiological system into a tractable form.

57

Compartment models are the most adopted approach for pharmacokinetics

58

modeling and allow the combination of dose and individual processes of absorp-

59

tion, distribution, metabolism and excretion (ADME) in a logical and straightfor-

60

ward manner (Rosenbaum, 2016). In this sense, the body is viewed as a system of

61

imaginary and homogeneous compartments, with no anatomical or physiologi-

62

cal representation and where drug concentration is uniform. Each compartment

63

correspond to groups of tissues with similar drug distribution rates, described by

64

mass balance equations.

65

In general, constant kinetic rates and other parameters are estimated using

66

curve fitting techniques from concentration-time profiles of drugs in the blood,

67

tissue, or other easily accessible body fluids such as urine, faeces or breast milk

68

(McMullin, 2005). The adjustment parameters obtained are mathematically rele-

69

vant, but they lack a physical meaning, which can limit the model as a predictive

70

tool, particularly for interspecies extrapolation. Another drawback of the model

71

is the difficulty, or even impossibility, to establish a straightforward relation-

72

ship between blood concentration-time profile and the intended site of action

73

concentration-time profile (Nestorov, 2003).

74

An alternative to the classical compartment models are the physiologically

75

based pharmacokinetics (PBPK) models. They also represent tissue or organs as 5

76

compartments, but describe the interconnected body compartments with perti-

77

nent physiological, biochemical and physicochemical information. Furthermore,

78

the route of drug administration respect their proper relationship with the over-

79

all physiology. Such features makes PBPK models a suitable tool for all types

80

of pharmacokinetic extrapolations, allowing the prediction of events in humans

81

from extrapolations of pre-clinical in vitro data or through intra-species extrapo-

82

lations from clinical data (Nestorov et al., 1998). A review of the PBPK modeling

83

approaches for drug research and analysis of advantages and limitations can be

84

found in Nestorov (2003).

85

The combination of dynamical modeling of cellular populations with pharma-

86

cological modeling has been used in a few studies. Among those are the works of

87

Scheiner et al. (2014) and Hambli et al. (2016) that incorporate a pharmacokinetic

88

(PK) model of the drug Denosumab, used in osteoporosis treatment, coupled to

89

the models of Pivonka et al. (2008) and Komarova et al. (2003), respectively, as a

90

pharmacodynamics model for quantifying the effect of denosumab on bone re-

91

modeling. Both works used a classical two compartment model, one for the sub-

92

cutaneous tissue and another for the blood serum. The results of Hambli et al.

93

(2016) showed an increase of BMD in the femoral neck, trochanter and total hip

94

compared with placebo in conformity with the results reported in the literature.

95

The model of Scheiner et al. (2014) allows reasonable simulations of osteoporotic

96

scenarios and pharmacological intervention through the administration of the

97

drug denosumab.

98

This research aims to integrate a physiologically-based pharmacokinetic model

99

(PBPK) with a pharmacodynamic model for bone remodeling based on the dy-

100

namics of bone cells populations. The modeling of cells populations describes 6

101

the signaling between osteoblasts and osteoclasts based on the model of Lemaire

102

et al. (2004), with the addition of the mechanical stimulus proposed by Scheiner

103

et al. (2013) and a pharmacological stimulus that will drive the bone remodeling

104

process. The PBPK model provides the drug concentration in the bone, which

105

allows the computation of a regulatory function that influences the dynamics

106

of bone cells population. The drug’s antiresorptive action is represented here

107

by a loss of efficacy of the active osteoclasts. We simulate this effect by adding

108

a regulatory function in the RANKL expression and increasing the osteoclasts

109

apoptosis rate. Therefore, the proposed model combines the advantages of phys-

110

iologically based pharmacokinetic models, whose compartments correspond to

111

exact anatomical or physiological entities, with a cellular population model that

112

better captures the complex osteoblast-osteoclast interactions in the microme-

113

chanical environment.

114

2. Materials and methods

115

This work proposes an unified model for bone remodeling, coupling the dy-

116

namics of bone cell populations and a pharmacokinetic model. In summary, the

117

procedure involves a finite element analysis performed in a three-dimensional

118

proximal femur. Geometry and material properties were modeled based on CT

119

images from previous work (Bahia et al., 2019). Once the structural analysis is

120

available, two stimuli to the bone cells are inputed to the remodeling process -

121

(i) a mechanical stimulus evaluated at the microscale, i.e. the microstrain energy

122

density; and (ii) a pharmacological stimulus, function of the drug concentration

123

on the bone. The strain energy density evaluated at the microscale is obtained

124

by a localization procedure (Mori-Tanaka scheme) performed in the nonhomo7

125

geneous bone medium. The concentration of the drug in the bone is provided by

126

a pharmacological model based on physiology, as can be seen in Fig. 2.

Figure 2: Overview of the bone remodeling procedure integrating the dynamics of bone cells populations and a physiologically based pharmacokinetics model. Mechanical loadings and local drug concentration are the stimuli in the bone remodeling process model. A micromechanical representation is used for the homogenization and scaling down of the strain energy densities for an extravascular bone matrix level. Art: Prepared by the author.

127

2.1. Physiologically based pharmacokinetic model

128

A physiologically based pharmacokinetic model describes the drug concen-

129

trations in the blood and other organs using physiological parameters such as

130

blood flow and organ volume as well as compartments representing different tis-

131

sues and organs. Here the PBPK model will provide the drug concentration at

132

every bone site, and is converted into a pharmacological stimuli in the dynamics 8

133

of bone cells populations model, represented by a system of ordinary differential

134

equations.

135

A schematic of the PBPK model structure is given in Fig. 3. As can be seen

136

in the model, the body is divided into a number of tissue compartments, each

137

one characterized by appropriate volume and flow rates. For simplicity, here

138

we adopt the following compartments: blood, gut, liver, kidney, rapid perfused

139

tissues, slowly perfused tissues and bone.

140

The temporal evolution of the drug concentration in each tissue/organ is de-

141

scribed with a mass balance equation that accounts for tissue blood flow rate

142

(Qi ) and the concentration difference between the arterial blood (cA ) and the ve-

143

nous blood (cVi ). Each tissue bone flow rate is assumed to be proportional to the

144

cardiac output (Q) and based on reference values for a 70 kg man (ICRP, 1975;

145

Yokley et al., 2006). Mean cardiac output is assumed to be 276 liters per hour.

146

Physiological values are adjusted to satisfy the requirement Q = ΣQi , where

147

i correspond to the tissues or organs: gut, kidney, liver, rapid perfused tissues,

148

slowly perfused tissues and bone. The volumes of each organ are considered

149

fractions of the body weight (BW) in [kg] as seen in Tab. 1.

150

The route of administration for a drug may be intravenous, or oral (gastroin-

151

testinal tract). In addition to these routes of administration, there are also pul-

152

monary (inhalation), subcutaneous (adhesives), muscular (injection) and oph-

153

talmic (eye drops). Drugs administered orally may be absorbed in different parts

154

of the gastrointestinal tract. Drugs absorbed through the oral cavity or rectum

155

bypass the liver. All other parts of the alimentary canal (especially the stom-

156

ach, duodenum, small and large intestines) carry blood (and drugs) to the liver

157

through the portal vein. There are also transporters in the mucosal cells of the gut 9

Figure 3: Physiologically based pharmacokinetic (PBPK) model. Qi terms are blood flow rates, c terms are concentrations. Arterial input is on the right and venous effluent in on the left. Source: Prepared by the author.

158

wall that pick up the drug from the gut and may excrete the drug back into the

159

intestine. Relatively high drug concentrations are processed by the liver during

160

the first passage of the drug. The liver can excrete the drug in the bile and then

161

into the duodenum. Therefore, the drug can be reabsorbed in an enterohepatic

162

cycle (Gieschke and Serafin, 2013).

163

The main organs for drug elimination are the liver and the kidneys. The liver

164

with its enzymatic equipment is the main metabolizing organ in the body. The

165

kidneys filter out wastes (including small drugs) from the blood and excrete them 10

166

with urine. Drug absorption by tissues depend on blood perfusion and tissue

167

permeability.

168

In the case of drugs used for treatment and prevention of osteoporosis, focus

169

of this investigation, they can be divided basically into two categories : anti-

170

resorptive agents or anabolic agents. Anti-resorptive agents reduce bone resorp-

171

tion, leading to increased BMD by several degrees. Estrogen, raloxifene, bispho-

172

sphonates and the human monoclonal antibody to RANKL are included in this

173

category. On the other hand, anabolic agents such as teriparatide (PTH1-34) and

174

total parathyroid hormone (PTH1-84), stimulate bone formation, thus increasing

175

BMD (Szulc and Bouxsein, 2011).

176

Most of the pharmacokinetic characteristics of bisphosphonates have been

177

derived from urinary excretion data since the plasma concentration following

178

therapeutic doses generally falls below the sensitivity limits of the quantification

179

assays (Porras et al., 1999). The bisphosphonate is eliminated from the plasma

180

by bone deposition and urinary excretion. For illustrative purposes, in rats it has

181

been found that 30 to 40% of 1mg/kg alendronate is eliminated in urine within 24

182

hours after dosing (Lin et al., 1992). About 60 to 70% of the dose is sequestrated to

183

the bone with short-term effect. In rats, the rate of renal excretion of alendronate

184

reached a maximum of about 25 mg/min/kg (Lin et al., 1991).

185

The mass balance equations for the PBPK model illustrated in Fig. (3) are

186

listed in Eqs. (1). The variables Vi , ci , Qi , t are volume, concentration and blood

187

flow in the tissue/organ i, and time, respectively. The indices i refers to blood, gut,

188

liver, kidney, rapidly perfused tissues (rp), slowly perfused tissues (sp) and bone,

189

respectively. Note that the variables Q, cV and cA represent the cardiac output

190

(blood flow), concentration in venous flow and concentration in arterial flow, 11

191

respectively. Indices A and V in each concentration expression refers to arterial

192

affluent and venous effluent, respectively, of an especific organ. For instance,

193

cA is the affluent concentration entering in each organ and cVbone is the effluent

194

concentration leaving the bone. Tab. (1) has the description and values of all the

195

symbols used in Eqs. 1a-1g.

dcblood = Q (cV − cA ) dt
dcgut Vgut = Qgut cA − cVgut dt dcliver Vliver = Qliver cA + Qgut cVgut − (Qliver + Qgut ) cVliver dt

Vblood

− CLhepato−bilis cVliver Fu
= Qkidney cA − cVkidney − CLrenal cVkidney Fu

dckidney dt
dcrp Vrp = Qrp cA − cVrp dt
dcsp = Qsp cA − cVsp Vsp dt dcbone = Qbone (cA − cVbone ) Vbone dt

Vkidney

(1a) (1b)

(1c) (1d) (1e) (1f) (1g)

196

The concentrations at t = 0 are cblood = cgut = cliver = ckidney = crp =

197

csp = cbone = 0. Furthermore, cA = cb , cVgut = cgut /Rgut , cVliver = cliver /Rliver ,

198

cVkidney = ckidney /Rkidney , cVrp = crp /Rrp , cVsp = csp /Rsp , e cVbone = cbone /Rbone .

199

Indices A and V in each concentration expressions refers to arterial and venous

200

effluent, respectively. Thus, the total venous concentrations is given by:

cV =

(Qliver + Qgut ) cVliver + Qkidney cVkidney + Qrp cVrp + Qsp cVsp + Qbone cVbone Q (2) 12

201

Finally the total blood flow is expressed by: Q = Qliver + Qgut + Qkidney + Qrp + Qsp + Qbone

(3)

202

The proposed PBPK model was fully written in MatLab 7.12 (R2011a) software

203

(MathWorks, Natick, Massachusetts, US) and run in a standard personal com-

204

puter (Intel Core i7 5500U CPU 2.4 GHz). The physiological parameters adopted

205

here combine data or are derived from different sources, such as O’Flaherty

206

(1991b); Pertinez et al. (2013) and are listed in Tab. (1).

207

2.1.1. Bone tissue modelling

208

The mineralized components of bone tissue are described by an adaptation

209

of the permeability rate limited tissue compartment model (Nestorov, 2003; Per-

210

tinez et al., 2013). According to this model, the bone is composed of a vascular

211

compartment, defined as the bone surface, and an extravascular compartment,

212

defined as the bone matrix. The bone remodelling process will govern the ex-

213

change between these two compartments. In fact, whereas bone is woven into

214

two main patterns (cortical and trabecular), each with different appearance and

215

characteristics, the bone tissue model presents a total of two compartments both

216

with two sub-compartments: cortical bone matrix, cortical bone surface, trabec-

217

ular bone matrix and trabecular bone surface. Expressing in mathematical terms,

218

the temporal evolution of the drug concentration in the vascularized bone sur-

219

face cvasi and in the extravascular bone matrix cbmi with index i assuming cortical

220

or trabecular types, can be modeled as:   dcvasi Qbone cA Qbone cvasi cvasi cbmi = − Rbone − kform + kres dt Vvasi Vvasi Vvasi Vbmi dcbmi cvasi cbmi = kform − kres dt Vvasi Vbmi 13

(4) (5)

Table 1: Physiological parameters used in the PBPK model.

Simbol

Units

Value

Description

Q

l/h

276

Cardiac output

Qgut

l/h

0.19Q

Blood flow to guta

Qkidney

l/h

0.19Q

Blood flow to kidneya

Qliver

l/h

0.06Q

Blood flow to livera

Qbone

l/h

0.19Q

Blood flow to bonea

Qrp

l/h

0.36Q

Blood flow to well perfused tissuesa

Qsp

l/h

0.125Q

Blood flow to slowly/poorly perfused tissues

BW

kg

70

Body weighte

Vgut

l

0.2895BW

Gut volumeb

Vkidney

l

0.004BW

Liver volumec

Vliver

l

0.025BW

Kidney volumec

Vrp

l

0.04BW

Well perfused tissues volumec

Vbone

l

0.026BW

Bone volumed

Vsp

l

4

Poorly perfused tissues volume

Rgut

n/a

0.9

Partition coefficient gut/plasmab

Rliver

n/a

0.47

Partition coefficient liver/plasmab

Rkidney

n/a

3

Partition coefficient kidney/plasmab

Rrp

n/a

0.9

Partition coefficient well perf. tissues/plasmab

Rsp

n/a

0.6

Partition coefficient poorly perf. tissues/plasmab

Rbone

n/a

2

Partition coefficient bone/plasma

CLrenal

n/a

0.03

Renal clearanceb

CLhepato bilis

n/a

0.03

Hepato-biliary clearanceb

Fu

n/a

0.73

Fraction unbound in blood b

a

Brown et al. (1997).

b

Pertinez et al. (2013).

c

Yokley et al. (2006).

d

O’Flaherty (1991b).

e

ICRP (1975).

14

221

where the indexes vas and bm refer to vascularized bone surface and extravas-

222

cular bone matrix, respectively. The terms kform and kres describe the bone for-

223

mation and bone resorption rates, respectively, and are both defined in the bone

224

cells populations model. Fig. 4 illustrates the bone tissue compartment model.

225

The steady-state values of total vascularized bone surface cvas and the total

226

extravascular bone matrix cbm will drive the dynamics of bone cells populations.

227

These total values are obtained by adding the cortical and trabecular contribu-

228

tions as follows:

229

cvas = cvascort + cvastrab

(6)

cbm = cbmcort + cbmtrab

(7)

where the indexes cort and trab refer to cortical and trabecular types, respectively.

Figure 4: Bone tissue compartment model. Source: Prepared by the author. 230 231

Bone physiological parameters such as the total bone volume, and the specific

232

surface and matrix volume fractions for the cortical and trabecular compartments 15

233

were extracted or derived from O’Flaherty (1991a) following the procedure of

234

Pertinez et al. (2013). The bone weight was calculated as 0.564 times total marrow

235

free dry bone weight considering a density of 1.92 g/ml for the marrow free dry

236

bone.

237

2.2. Cells populations model

238

Bone remodeling is a biological process of bone adaptation that involves se-

239

quential osteoblast-mediated bone formation and osteoclast-mediated bone re-

240

sorption at the same location, acting to maintain bone homeostasis, repairing

241

small defects and replacing bone matrix (Burr and Allen, 2013). A complex net-

242

work regulates the interactions between bone cells. The signaling pathway RANK-

243

RANKL-OPG has been shown to be a major regulator of bone remodeling act-

244

ing on osteoclastogenesis (Hofbauer et al., 2004; Martin, 2004). RANKL is a cy-

245

tokine (membrane bound factor) produced by osteoblasts and stromal cells that

246

binds to its receptor RANK expressed by osteoclasts, and induces both differ-

247

entiation and activation of osteoclasts. Osteoprotegerin (OPG), also produced

248

by stromal/osteoblastic cells, binds to RANKL to prevent RANKL stimulation by

249

RANK, inhibiting bone resorption (Pivonka et al., 2008).

250

Imbalances in the RANK-RANKL-OPG pathway are closely associated with

251

diseases such as osteoporosis, osteopenia, and tumors, among others. It is spec-

252

ulated that an overproduction of RANKL is associated with the development

253

of postmenopausal osteoporosis (Miller, 2009). Furthermore, the transforming

254

growth factor (TGF-β) stimulates osteoblastic differentiation and promotes bone

255

formation and hormones like the parathyroid hormone (PTH) acts by inhibiting

256

OPG and enhancing RANKL expression.

257

In this context, mathematical modeling provides a powerful tool to study the 16

258

interactions between bone cells, their biochemical regulators, allowing the un-

259

derstanding of the bone remodeling process in distinct scenarios. Lemaire et al.

260

(2004) proposes a mathematical model for the temporal evolution of the bone

261

cell populations through a system of differential equations. The model includes

262

interesting features like the osteoblast-osteoclast coupling, the RANK-RANKL-

263

OPG pathway, and the antagonistic effects of PTH administration on the control

264

of OPG and RANKL expression by active osteoblasts. These features enabled

265

the assessment of certain bone metabolic dysfunctions. In addition, the model

266

includes the effect of growth factor TGF-β on the differentiation of precursor

267

osteoblasts.

268

Here, some modifications were applied to the original model by Lemaire et al.

269

(2004) to incorporate a mechanical stimulus based on strain energy density eval-

270

uated at the microscale (Scheiner et al., 2013) and the pharmacological regulation

271

of an anti-resorptive drug on the bone matrix (Ross et al., 2017). The dynamics

272

of the proposed bone cells populations obtained, built upon the work of Lemaire

273

et al. (2004) are defined in Eqs. 8-18. The variables OBp , OBa and OCa refer to

274

pre-osteoblasts, active osteoblasts and active osteoclasts, respectively. BA and

275

BI are the antiresorptive drug concentrations in the bone matrix and plasma, re-

276

spectively. The behavior of bone cells is regulated by the binding of cytokines to

277

the receptors on the cells surface. The chemical kinetics of reactions can be rep-

278

resented by the law of mass action. A graphical scheme of the model is presented

279

in Fig. 5.

17

Figure 5: Sketch of the bone cells populations model. Microscopic strain energy density (microSED) and drug concentration plays a regulatory role in the dynamics of cells populations. These stimuli are added to the previous work of Lemaire et al. (2004). Artwork: Miguel T. Bahia.

dOBp dt dOBa dt dOCa dt dBA dt

= DR πC + POBp OBp Πε − =

DB OBp πC

DB OBp − kB OBa πC

(8) (9)

= DC πL − (1 + k7 BA)DA πC OCa

(10)

= k8 cdrug (1 − BA) − k9 BA

(11)

280

where cdrug is the local drug concentration that reaches the bone obtained from

281

the PBPK model. For simplicity, we adopted the steady-state value of the total

282

extravascular bone matrix, i.e., cdrug = cbm . BA is the drug amount that binds

283

to the bone matrix. DB is a factor of proportionality and πC is the influence of

284

TGF-β, expressed by: DB = f0 dB πC = 18

C + f0 C S C + CS

(12)

285

RANKL induces the differentiation of active osteoclasts which release TGF-β

286

causing osteoclastic apoptosis and controlling osteoblast differentiation.

287

The drug concentration in the bone obtained from the PBPK model is used

288

here to trigger an antiresorptive effect through a pharmacoregulation function

289

ΠBA .

290

291

The RANKL influence, πL is obtained by:   IL + PRLεbm k3 KLP .πP .OBa  P  1 + ΠBA + πL = KO k4 1 + k3 + k1 rL k4 k2 k0 πP +I0 where πP is the fraction of occupied PTH receptors which is expressed by: πP =

292

293

(13)

IP kP IP kP

+ +

SP kP k6 k5

(14)

and the term PRLεbm corresponds to the variation of RANKL related to changes in strain energy. It is activated only when in disuse and is expressed by: ! wεbm PRLεbm = κ 1 − wεbmst

(15)

294

where the adjustment constant κ is null for wεbm ≥ wεbmst . and positive other-

295

wise. The incorporation of this non-null term mathematically expresses that any

296

deformation would limit the external production of RANKL by osteocytes and

297

cap marrow stromal cells, and would be totally stopped if the critical value of

298

SED is reached. Physiologically, the mechanotransductive activity of osteocytes

299

located in lacunar pores of the extravascular bone matrix would be signaled in at

300

least three ways: direct deformation of osteocytes, hydrostatic pressure or even

301

the presence of fluid flow. As Scheiner et al. (2013) points out, all of these sit-

302

uations are associated with a decrease in the RANKL/ OPG rate expressed by

303

osteocytes. 19

304

The term ΠBA is a drug-regulatory function that simulates the effect of the

305

external administration of an antiresorptive drug on the population of cells. The

306

term was added to the original expression proposed by Lemaire et al. (2004) for

307

the RANKL concentration equation, being expressed by: ΠBA = kBA BA

308

(16)

where BA is given by Eq. 11 and kBF is an adjustment constant.

309

The population of preosteoblasts OBp increases with the differentiation of

310

uncommitted progenitors osteoblast cells with a maximal differentiation rate DR ,

311

promoted by the growth factor TGF-β.

312

Uncommitted progenitor osteoblasts represent the initial stage of maturation

313

of the osteoblast lineage, represented by a constant pool of mesenchymal stem

314

cells located in the bone marrow, spleen and liver (Gowder, 2012). As pointed

315

out by Lemaire et al. (2004), this category corresponds to a composite of several

316

phenotypes that share similar characteristics.

317

On the other hand, the population of precursor osteoblasts decreases due

318

to the differentiation of the precursor osteoblasts into active osteoblasts which

319

occurs at a maximum rate.

320

Following the proposal of Scheiner et al. (2013), the mechanical stimulus is

321

modeled adding a term in Eq. 8 that regulates the proliferation of osteoblasts

322

as a function of the strain energy density evaluated locally. The function Πε

323

is the mechanoregulation and PR is the constant indicating the proliferation of

324

osteoblasts according to the strain energy density:

Πε = Πεst

1+λ 20

!! wεbm −1 wεbmst

(17)

325

where Πεst is the equilibrium value of Πε , wεbm is the value of local strain en-

326

ergy density and wεbmst is the equilibrium value. Parameter λ is an adjustment

327

constant that is null to wεbm ≤ wεbmst and positive otherwise. In the presence of

328

bone diseases such as osteoporosis or even during a treatment with antiresorp-

329

tive agents this parameter should be recalibrated since it significantly affects the

330

bone response to mechanical stimuli. Current studies suggest that precursor os-

331

teoblasts may respond to mechanical stimuli, either matrix deformations and/or

332

fluid flow, and osteocyte-derived biochemical signaling (i.e. sclerostin). In the

333

present work, a phenomenological activating function was adopted, following

334

the proposal of Scheiner et al. (2013), as seen in Eq. 17. As shown in Eqs. ? and

335

?, the maximum proliferation rate of precursor osteoblasts is reached when the

336

mechanoregulator term equals to 1, which would provide sufficient mechanical

337

activation of osteoblasts. According to the studies of Jones et al. (1991) and Kas-

338

par et al. (2002), low straining levels would reduce the proliferation rate by 25 and

339

50%, which is simulated by the adoption of a minimum value of the mechanoreg-

340

ulatory function Πεst related to the limit value wεbm . Above this value the in-

341

creased proliferation of precursor osteoblasts is activated. Thus, the function

342

provides intermediate values between the maximum and minimum proliferation

343

rates.

344

The model’s dead zone is defined in the mechanical and pharmacological

345

stimuli expressions. In the case of mechanical stimulation, Eq. 17 provides null

346

values when the strain energy density is within the corresponding range. When

347

the drug concentration is zero, the expression of the pharmacological stimulus

348

produces null values.

349

The active osteoblast population OBa increases according to the pre-osteoblast 21

Table 2: Parameters used in the simulation of the bone remodeling model. Values are extracted from Lemaire et al. (2004), except for ki (i = 7, ..., 12) and kBA , that were adjusted by the authors Symbol

Units

Value

Description

5 × 10−3

Number of OC to get half differentiation flux

DA

day

−1

0.7

Rate of osteoclast apoptosis caused by TGF−β

dB

day−1

0.7

Cs

pM

Differentiation rate of responding osteoblasts

DC

pM.day

−1

2.1 × 10

Differentiation rate of osteoclast precursors

DR

pM.day−1

7 × 10−4

Differentiation rate of osteoblast progenitors

−3

f0

n/a

0.05

Fixed proportion

IL

pM.day−1

0 − 106

Rate of administration of RANKL

IO

pM.day−1

0 − 106

Rate of administration of OPG

IP

pM.day

6

0 − 10

Rate of administration of PTH

10

Fixed concentration of RANK

10−2

Rate of OPG-RANKL binding

−1

K

pM

k1

pM−1 .day−1

10

Rate of OPG-RANKL unbinding

k3

pM−1 .day−1

5.8 × 10−4

Rate of RANK-RANKL binding

k4

pM.day−1

1.7 × 10−2

Rate of RANK-RANKL unbinding

k2

day

−1

0.02

Rate of PTH binding with its receptor

k6

day−1

3

Rate of PTH unbinding

kB

day−1

0.189

KLP

pM/pM

kO

day−1

KOP

day

2 × 10

Minimal rate of production of OPG per cell

kP

day−1

86

Rate of elimination of PTH

rL

pM.day−1

103

Rate of RANKL production and elimination

SP

pM.day

k5

pM .day −1

−1

−1

Rate of elimination of active osteoblasts

3 × 10

Maximum RANKL attached on each cell surface

0.35

Rate of elimination of OPG

6

5

250

Rate of synthesis of systemic PTH

λ

n/a

1.2

Parameter of anabolic adjustment

κ

n/a

5 × 102

Parameter of RANKL inibition

Πεst

n/a

0.5

Equilibrium of mechanoregulation function

wεbmst

Pa

0.78

SED inferior limit

k7

n/a

0.3

−1

Hill function parameter for apoptosis of OCa

6 × 10

Bisphosphonate transfer from plasma to bone

k9

1/h

5 × 10−4

Bisphosphonate on bone clearance rate

k10

1/h

25 × 10−5

First order BMD rate OCa /OBa

k11

1/h

1.5 × 10

BMD balance

k12

n/a

0, 7979

Baseline value of the OCa /OBa

kBA

n/a

1

Parameter of the Hill function for drug regulation

k8

pM.day

−1

−3

−4

22

350

differentiation, which occurs with the maximum differentiation rate DB and is

351

inhibited by the TGF-β. The population is reduced by the apoptosis of active

352

osteoblasts, which follows the apoptosis rate kB .

353

Antiresorptive drugs such as bisphosphonates bind to the hydroxyapatite in

354

the bone and inhibit resorption by osteoclasts. In the model, the loss of func-

355

tion of osteoclasts is simulated by a change in RANKL concentration in Eq. 13

356

and by an increase in the rate of apoptosis of active osteoclasts in Eq. 10. The

357

bisphosphonate attaches to the bone in a second-order forward, first-order re-

358

verse reaction with hydroxyapatite attachment sites (Ross et al., 2017), resulting

359

in values in the range [0, 1]. In Eq. 11, the drug concentration in the bone ma-

360

trix BA follows Ross’s proposal. The plasmatic anti-resorptive agent BI attaches

361

to hydroxyapatite cristals in bone with a plasma-to-bone transfer rate k8 and is

362

eliminated at a rate k9 .

363

The model proposed by Ross et al. (2017) incorporates a reference biomarker

364

for the bone mineral density (BMD) into the system of differential equations (Eqs.

365

8-11), here denoted as the dimensionless variable BM in the Eq. 18. The evolution

366

of BMD is governed by the difference between the osteoblast/osteoclast ratio and

367

its reference value, and by the difference between the biomarker and a baseline

368

value: dBM = k10 dt



 OBa − k12 − k11 (BM − 1) OCa

(18)

369

The BM biomarker present in Eq. 18 allows the analysis of the bone mass in re-

370

lation to a reference value stipulated by the user. In this equation, 1 was adopted

371

as a reference, but other values can be adopted depending on the clinical con-

372

text. Note that, in the absence of treatment, the biomarker for BMD returns to

373

the baseline value 1. The inclusion of Eq. 18 in the model makes it more flexible, 23

374

enabling site-specific BMD changes, according to the clinical data (Ross et al.,

375

2017).

376

The model parameters are summarized in Table 2. A more detailed explana-

377

tion of each parameter is found in Lemaire et al. (2004); Scheiner et al. (2013);

378

Pivonka et al. (2008).

379

2.3. Microscale strain energy density

380

The strain energy density evaluated at the microscale (microSED) is adopted

381

as a mechanoregulatory stimulus for osteocytes to drive the bone remodeling

382

process (Jacobs, 1994; R¨uberg, 2004). According to the magnitude of this energy,

383

the algorithm distinguishes remodeling zones with different behaviors.

384

The microscopic evaluation of strain energy densities follows the proposal

385

of Scheiner et al. (2013), which considers the inhomogeneous constitution of

386

the bone and the resulting inhomogeneity of its mechanical properties (Zaoui,

387

2002). The transition between macroscale and microscale involves techniques

388

of homogenization that originates from the continuum micromechanics. The

389

bone presents a hierarchical structure quite complex relative to its microscopic

390

arrangements and mechanical properties. These procedures represent a viable al-

391

ternative to obtaining effective mechanical properties of the material to be used

392

in the structural analysis. In this direction, one can cite the work of Fritsch and

393

Hellmich (2007) that proposes a multiscale material model for the bone consid-

394

ering its hierarchical structure in a multi-step homogenization procedure, re-

395

sulting in the estimation of the elastic properties of the bone at the nano and

396

sub-microscale.

397

For the sake of simplicity, we chose a characteristic length for representative

398

volume element (RVE) that allows the bone to be reasonably represented as a 24

399

biphasic composite of extravascular bone matrix with vascular pore space in be-

400

tween, as shown in Fig. 6. Each constituent has a corresponding volume fraction

401

fi , defined as fi =

402

the total volume, obeying the relation:

Vi , Vtotal

with Vi as the volume of the constituent i and Vtotal as

fbm + fvas = 1

(19)

403

with fbm as the extravascular bone matrix and fvas as vascular pores space. Due

404

to the inhomogeneities, the volume fractions of each of these phases vary accord-

405

ing to the type and bone site. In the case of cortical bone, the volume fraction

406

of the extravascular bone matrix lies within the range 0.7 - 0.95, whereas for the

407

trabecular bone the values are between 0.1 and 0.5 (Scheiner et al., 2016; Cooper

408

et al., 2007; Padilla et al., 2008; Boutroy et al., 2005, 2011).

409

410

It is possible to express fvas (or fbm ) in terms of the apparent density ρ = mtotal , Vtotal

with mtotal as total mass. Thus:

fvas = 1 −

Vbm (mtotal /ρ)

(20)

411

Assigning zero mass to the vascular pore space, i.e., mtotal = mbm , the expression

412

is reduced to: fvas = 1 −

ρ ρ or fbm = ρt ρt

(21)

413

where ρt = mbm /Vbm is the density of an imaginary bone specimen without

414

vascular pores (R¨uberg, 2004).

415

The temporal evolution of the bone constituent volume fraction is incorpo-

416

rated into the bone cells populations model. Thus, the volume fraction of ex-

417

travascular bone matrix decreases due to bone resorption and increases due to 25

418

the bone formation. In mathematical terms, this can be represented by the equa-

419

tion: df bm = kform OBa − kres OCa dt

(22)

420

with kform as the formation rate (bone mass formed per unit of time) and kres as

421

the resorption rate (bone mass resorbed per unit of time). Conversely, the volume

422

fraction of vascular pore space decreases due to bone resorption and increases

423

due the bone formation, which results in the condition: dfvas dfbm =− = −kform OBa + kres OCa dt dt

(23)

424

Both formation and resorption rates, kform and kres respectively, are consid-

425

ered species-specific cellular properties, and although known to be dependent on

426

health conditions, here they are considered invariant in time. kform is calibrated

427

so that the previously defined bone resorption rate and cells populations values

428

for stationary states imply in a balanced bone turnover (in the EDO system, we

429

have dOBp /dt = dOBa /dt = dOCa /dt = 0) thus considering the bone matrix

430

volume fraction fbm constant in Eq. 23, kf orm = kres OCa /OBa .

431

432

From the perspective of the continuum micromechanics, the homogenized fourth order stiffness of the bone Chom bone , is: Chom bone =

X

fr Cr : Aest r

(24)

r 433

where Cr is the microscopic fourth-order stiffness tensor of the constituent r,

434

with r = vas, bm, and Aest r is the estimation of the corresponding fourth-order

435

strain concentration tensor, relating the macroscopic and microscopic second or26

Figure 6: Micromechanical representation of the bone as a two-phase composite material of ”vascular” pores and ”extravascular” bone matrix with volume fractions fvas and fbm , respectively. Vascular pores host osteoblasts and osteoclasts at various differentiation stages and a large variety of biochemical factors. Artwork: Miguel T. Bahia.

436

der strain tensors. Aest r can be estimated by solving the classical matrix-inclusion

437

problem of Eshelby (Eshelby, 1957), by means of the Mori-Tanaka scheme:

Aest = r

 −1  I + Pbm : r : (Cr − Cbm ) )−1 ( X   −1 fs I + Pbm s : (Cr − Cbm )

(25)

s 438

where the index r denotes the phase, and the sum over the index s include both

439

phases, s = vas, bm. In addition, I is the fourth-order unit tensor with the com-

440

ponents defined by the Kronecker delta, δij = 1 for i = j and zero otherwise,

441

being Iijkl = 1/2 (δik δjl + δil δjk ), and Pbm r is the Hill tensor of phase r inserted

442

in the matrix with rigidity Cbm .

443

Volume fractions are obtained from the temporal evolution of the cellular 27

444

population model. Stiffness tensors are extracted from the work of Scheiner et al.

445

(2013).

446

A linear elastic constitutive law relates the macroscopic stress tensor acting

447

on the bone and the corresponding macroscopic strain tensor, via the macro-

448

scopic stiffness tensor (Eq. 24), i.e.: hom −1 : Σbone Σbone = Chom bone : Ebone ⇐⇒ Ebone = (Cbone )

(26)

449

The strain tensor of the extravascular bone matrix, corresponding to the

450

macroscopic loading Σbone , can be expressed with the stress concentration tensor

451

Aest bm defined in Eq. 25, that is: est hom −1 εbm = Aest : Ebone ] bm : Ebone = Abm : [(Cbone )

452

453

(27)

The microscopic strain energy density in the extravascular bone matrix is expressed by: 1 Ψbm = εbm : Cbm : εbm 2

(28)

454

The microscopic SED Ψbm obtained by the homogenization technique may

455

now be used as a mechanoregulatory function providing a representative scalar

456

value at the microscale associated with the three-dimensional macroscopic stress

457

state.

458

The bone cells populations model was implemented incorporating a FE frame-

459

work for the three-dimensional case. Each interaction of the structural analysis

460

calculates the local micro strain energy density for each node of the discretized

461

mesh, thus solving the system of differential equations that governs the tempo-

462

ral evolution of the bone cells populations. A new population of bone cells is 28

463

calculated at each node and the volume fractions are updated at each iteration.

464

The temporal evolution of the bone matrix volume fraction (Eq. 22) provides the

465

temporal evolution of apparent bone density ρ from Eq. 21. The homogenized

466

stiffness tensor Chom bone is updated for each node of the finite element mesh using

467

the new volume fractions. Note that in the first iteration, a phenomenological

468

law (Eqs. 29-30) (described in Subsection 2.4) relating the elasticity modulus and

469

the ash densities from the CT images is used in the structural analysis. The next

470

iterations use the updated homogenized stiffness tensor in the structural analy-

471

sis.

472

Bone is largely complex, both physiological and morphologically. Depending

473

on the observation scale (from macroscopic to nanoscopic), different constituents

474

with distinct mechanical properties are observed: hydroxyapatite crystals, col-

475

lagenous proteins, water, among others. Thus, to obtain the effective properties

476

of the composite (particularly the stiffness matrix) the Mori-Tanaka homogeniza-

477

tion technique was used. The actual modulus of elasticity of the bone does not

478

depend solely on ash density as the phenomenological relationship may suggest,

479

but it was adopted here for the sake of simplicity. Literature offers several ex-

480

perimental relationships between density and elastic modulus for different bones

481

(tibia, femur, spine) and different regions of the bone. Some examples are Carter

482

and Hayes (1977); Morgan et al. (2003); Snyder and Schneider (1991); Kaneko et al.

483

(2004); Keller (1994). A review of these relationships can be found in Helgason

484

et al. (2008).

485

2.4. Finite Elements Model

486

In order to apply the proposed procedure, a FE model was built from CT

487

images of a femur. The initial densities were derived from the CT images and 29

488

inputed to the FE model, as shown in Fig. 17a. We adopted a linear relationship

489

between pixel values [0.255] and bone densities.

490

There are several relationships between Young’s modulus and apparent bone

491

densities in the literature (Carter and Hayes, 1977; Morgan et al., 2003; Snyder

492

and Schneider, 1991; Kaneko et al., 2004; Keller, 1994). A review can be found in

493

Helgason et al. (2008).

494

As observed previously, in the first iteration the elasticity modulus is given

495

as a function of the ash densities derived from CT images using the relation

496

proposed by Keller (1994) and validated by Yosibash et al. (2007) and Trabelsi

497

et al. (2009):

498

Ecort = 10200ρ2.01 ash [M P a]

(29)

Etrab = 5307ρash + 469 [M P a]

(30)

Trabecular and cortical bone types are distinguished with the condition ρash >

499

0.6g/cm3 for cortical bone and ρash ≤ 0.6g/cm3 for trabecular bone (Trabelsi

500

et al., 2014).

501

The model was fully implemented in Matlab and post-processed in Paraview

502

5.0.0. The procedure for reconstruction and mapping of non-homogeneous isotropic

503

material properties is described in more details in Bahia et al. (2016).

504

The boundary conditions were based on Carter et al. (1987) and represent the

505

joint reactions and abductor muscle forces for a standard gait cycle. For each

506

cycle of FE analysis, three loading cases corresponding to the maximal adduc-

507

tion limits, maximum abduction and average support of the gait cycle were an-

508

alyzed sequentially (Fig. 8). The orientations and magnitudes of the forces were 30

(a) 3D femur mesh.

(b) Mapping of densities.

Figure 7: Discretization of the reconstructed femur in a finite element mesh with 3846 nodes and 16689 tetrahedral elements (Fig. 17a) and assignment of gray values to mesh elements (Fig. 17b). Source: Prepared by the author. Table 3: Loads used in the FE model referring to a standard gait. Orientation angles to the frontal and sagittal planes (Hambli et al., 2016; Carter et al., 1987).

Load

Orientation [deg]

Orientation [deg]

Frontal

Frontal

Sagittal

Sagittal

Case 1

2317

24

6

703

28

15

Case 2

1158

−15

35

351

−8

9

Case 3

1548

56

−20

468

35

16

509

adapted from Hambli et al. (2016) and are listed in Tab. 3. The concentrated loads

510

proposed by Hambli et al. (2016) were transformed into equivalent proportional

511

distributed loads on the femoral head and greater trochanter.

512

The strain energy density (microSED) from the combined loading cases was

513

defined as the highest nodal value. In other words, considering the three finite 31

(a) Load case 1

(b) Load case 2

(c) Load case 3

Figure 8: Load cases adopted in the 3D FE analysis. Each load specified in Tab. 3 was proportionally distributed onto the femur head and greater trochanter to simulate the stress on joint during gait. Source: Prepared by the authors.

514

element structural analyses performed for each load case, we adopted the critical

515

value (maximum absolute value) of equivalent stress obtained from the three

516

simulations in each mesh node in the bone remodeling algorithm. This choice

517

produced the maximum mechanical stimulus from the combination of the three

518

loading cases. This is justified because the static analysis does not incorporate

519

the loading history and load cycles, which usually arise in real daily activities.

520

521

The pharmacodynamic (PD) model coupling pharmacokinetics and bone remodeling proposed could be summarized in Fig. 9.

32

Figure 9: Diagram for the proposed unified bone remodeling model. A pharmacokinetic model provides the plasma concentration of the drug that contributes to the evolution of bone cells populations. Source: Prepared by the author.

522

3. RESULTS

523

The proposed PBPK model for bone seeking agents was solved by the 4th or-

524

der Runge-Kutta method and implemented in Matlab. The dosing regimen cho-

525

sen was 10 mg/day intravenous infusion of antiresorptive drug over 4 hours (20

526

h intervals). In clinical practice, intravenous infusion offers an alternative ther-

527

apeutic option with less gastrointestinal disturbances than orally administered

528

bisphosphonates (Cremers et al., 2002). 33

529

Figs. 10 and 11 illustrate the time course of drug concentration in each tissue

530

or organ for the 10 mg regimen by intravenous infusion for 4 hours per day. Note

531

that each dose (10 mg) is divided by infusion time (4 h). The organs included in the

532

model are blood, liver, kidney, stomach, well-perfused tissues (rapid perfusion),

533

poorly perfused tissue (slow perfusion), and bone. Fig. 11 presents the temporal

534

evolution of the drug separately for each organ.

535

The intravenous infusion is modeled by periodically adding the drug dose,

536

mathematically expressed here by a step or Heaviside function. The partition

537

coefficient of the kidney is greater than the coefficient of the other organs, so the

538

concentration of the drug reaches higher peaks in this organ, as can be observed

539

in Fig. 11.

540

The bone compartment includes the vascular bone surface and the extravas-

541

cular bone matrix for both cortical and trabecular bone. For this regimen, the

542

drug is not completely eliminated from the bone. In this case, the drug con-

543

centration is divided in more compartments as presented in Sec. 2.1.1. In the

544

extravascular bone matrix, a portion of the drug is not completely eliminated

545

before the next infusion, and accumulates in the tissue.

546

Fig. 12 shows the drug distribution in the four bone compartments: the cor-

547

tical and trabecular bone and their respective subdivisions, the vascular bone

548

surface and the extravascular bone matrix. The model makes it possible to cali-

549

brate the drug distribution by modifying parameters such as bone partition co-

550

efficient, drug bioavailability and rates of bone formation and resorption, given

551

by coefficients kform and kres , related to the process of bone remodeling.

552

The physiological parameters used in the simulation are listed in Tab. 1. The

553

drug concentration in the bone is given as micrograms of drug per gram of bone 34

Figure 10: Implementation of the proposed PBPK model. A dosage regimen of 10 mg of drug / day by intravenous infusion over a period of 4 hours. The chart shows the rate of infusion and amounts of drugs in each organ / tissue. Source: Prepared by the author.

554

ash and is defined as cdrug = 0.01 mcg drug/g bone ash. This value corresponds to

555

the steady state drug concentration obtained in the trabecular bone matrix with

556

the PBPK model for the 10 mg daily drug regimen intravenously for 4 h. Note

557

that the most of the physiological parameters of the PBPK model used here were

558

extracted from the literature. However, due to the lack of specific experimental

559

data, some parameters, such as bone formation and resorption rates and bone

560

partition coefficients, needed to be arbitrated by the author. It is worth noting

35

Figure 11: Concentration of drug in each organ / tissue for a dosage regimen of 10 mg of drug per day, intravenously for 4 hours.

561

that the value adopted for the drug concentration in the bone has the purpose of

562

illustrating the procedure, but needs to be validated experimentally.

563

The bone cells population model represented by Eqs. 8-18 is ilustrated on Figs.

564

13 and 14. The model was solved by numerical integration 4th order Runge-Kutta

565

method. In the first simulation (Fig. 13), the strain energy density was defined

566

as 0.13 Pa and there was no bisphosphonate stimulation. Around the 180th day

567

of iteration the population of preosteoblasts, active osteoblasts and osteoclasts

568

reached the steady state condition. At the end of 200 days the BMD biomarker 36

Figure 12: Drug distribution in the trabecular and cortical bones for a dosage regimen of 10 mg of drug /day intravenously for 4 hours.

569

had an increase of 1.98%.

570

Note that cell populations maintain the steady state value when disturbance

571

is removed from the system of linear differential equations. To analyze the im-

572

pact of the system perturbations such as the inclusion of RANKL or OPG, it is

573

more significant to evaluate the system’s behavior in terms of percent change

574

from the current state relative to the initial reference values. After the stimu-

575

lus, as the pharmacological and/or mechanical, a new population configuration

576

is obtained, and the disturbed system steady state value will be different from

577

the initial steady state value, as shown in Fig. 13. Removing the system dis-

578

turbance, the cells populations return to the initial values, as shown by Lemaire

579

et al. (2004).

580

The second simulation aims to assess the impact of the anti-resortive agent 37

Figure 13: Dynamics of bone cell populations formed by pre-osteoblasts, active osteoblasts and active osteoclasts, and bone mass density (BMD) evolution. Mechanical stimulation by SED. Without anti-resorbtive agent.

581

on bone cells population and BMD biomarker, as shown in Fig. 14. The me-

582

chanical stimulus was maintained. The antiresorptive effect was simulated as

583

an increase in the rate of apoptosis of active osteoclast and the inclusion of a

584

regulatory function in RANKL expression. We observed that the increase of os-

585

teoclast apoptosis rate was not sufficient to increase the rate of differentiation

586

of precursor osteoblasts in the model, so we added another term (a regulatory

587

function given by Eq. 16) to the RANKL expression to produce the desired effect.

588

Repeating the simulation with the pharmacological stimulus after the 20th day,

589

we observed an increase of 3.28% in BMD biomarkers, slightly higher than that

38

590

observed in the previous simulation. The sensitivity of the BMD biomarker to the

591

mechanical stimulus, that changes the rate of differentiation of pre osteoblasts

592

was higher than its sensitivity to the pharmacological stimulus, that acts on the

593

apoptosis of osteoclasts and the expression of RANKL in the model. This suggests

594

that other pathways could be considered in the expression of pharmacological

595

intervention, including adjustments in mechanoregulation function. It is reason-

596

able to assume changes in sensitivity to mechanical stimulus in the presence of

597

certain bone diseases such as osteoporosis.

Figure 14: Dynamics of bone cell populations formed by pre-osteoblasts, active osteoblasts and active osteoclasts, and bone mass density (BMD) evolution. Mechanical stimulation by SED. With anti-resorptive agent.

598

The expression that defines the differentiation of pre-osteoblasts is a function 39

599

of the concentration of active osteoclasts. Thus, the increase in active osteoclasts

600

apoptosis does not increase the population of active osteoblasts, but reduces

601

the rate of differentiation of precursor osteoblasts. The increase in osteoclast-

602

osteoblast ratio was observed when the drug-regulatory function (Eq. 16) was

603

included in Eq. 13.

604

3.1. Postmenopausal osteoporosis treatment simulation

605

Postmenopausal osteoporosis alters bone microarchitecture, resulting in in-

606

creased porosity, reduced mechanical resistance and increased risk of fractures.

607

In order to verify the validity of the model in this scenario, we have simulated

608

postmenopausal osteoporosis and the effect of a drug treatment. The increased

609

porosity can be simulated in different ways in our proposed model. For instance,

610

by reducing the minimum rate of OPG production per cell, or modifying the rate

611

of administration of PTH and/or RANKL. For simplicity, here we chose to reduce

612

the minimum rate of OPG production per cell, according to Lemaire et al. (2004).

613

Considering that the bone adaptive response to mechanical stimuli is drasti-

614

cally affected by the changes in the RANK-RANKL-OPG system caused by osteo-

615

porosis, the anabolic factor λ in the mechanoregulation equation was calibrated

616

to favor osteoclast activity while decreasing osteoblast activity, which results in

617

greater bone resorption. Changes in anabolic factor follow the work of Scheiner

618

et al. (2013).

619

Bisphosphonates acts on osteoporosis by preventing the fixation of osteo-

620

clasts in the bone matrix and increasing osteoclastic apoptosis, thus inhibiting

621

bone resoption.

622

Fig. 15 illustrates the changes in the bone matrix volume fraction over time,

623

with and without antiresorptive agents. Pharmacological treatment was initiated 40

Figure 15: Temporal evolution of the bone matrix volume fraction with and without antiresorptive treatment. Pharmacological treatment started after the 200◦ day. Initial volume fraction fbm equal to 0.95.

624

after the 200◦ day. The minimum rate of OPG production per cell, represented by

625

the term koP , was reduced by half, i.e. 1.105 pM.day−1 /pM cells, and the anabolic

626

factor λ was set to 1.2. Pharmacological treatment was simulated by assigning

627

non-zero values to the variable cdrug that corresponds to the anti-resorptive drug

628

concentration in the bone, Eq. 11. The temporal evolution of bone matrix volume

629

fractions fbm for cdrug = 0 (no pharmacological stimulus), cdrug = 0.01, 0.1 and

630

1.0 mg/g bone ash is given in Fig. 15. The initial value for bone matrix volume

631

fraction fbm was 0.95. The increase in the pharmacological stimulus produced

41

632

a smaller decrease of the volume fraction, which is expected for this type of

633

treatment. The simulation was investigated for a period of 1200 days.

634

When the cell concentration reaches a plateau, the volume fraction starts de-

635

caying linearly. This occurs because the rates of bone formation and resorption,

636

kres and kform , which control the temporal variation of the bone matrix volume

637

fraction were considered constant in this work.

638

3.2. Bone remodeling simulation

639

In this section, we present the implementation of the proposed procedure as

640

shown in Fig. 9. Due to practical and ethical reasons, obtaining in vivo informa-

641

tion about drug disposition in human bones is challenging. Specimen collection

642

can damage bone tissue and is rarely conducted in the clinical setting Stepensky

643

et al. (2003). Drug disposition also depends on the type of anti-resorptive agent,

644

with significant differences even among bisphosphonates, such as alendronate

645

and etyledronate Masarachia et al. (1996). For this reason, the authors adopted

646

only one reference value for the simulations, in order to qualitatively evaluate

647

the procedure. Drug concentration in the bone used as the pharmacological stim-

648

ulus was cdrug = 0.5mcg/g bone ash. This value roughly corresponds to the drug

649

concentration in the bone matrix obtained for the 10 mg daily dosing regimen

650

via 4 h intravenous infusion as shown in Fig. 11.

651

Each procedure iteration involves a structural analysis using finite element

652

method, numerical integration of the system of differential equations that gov-

653

erns the dynamics of cells populations, and update of the mechanical properties.

654

The temporal evolution of the cells populations is performed in cycles of 120

655

days. Figure 17 shows the distribution of densities in a longitudinal section of

656

the reconstructed femur for the initial and final configurations after 4 iterations 42

(a) Longitudinal section plane.

(b) Cross section plane.

Figure 16: Cutting planes adopted in the visualization. The length of femur piece is about 40 cm and the cross-section was made at a height of 10 cm measured from the lower end.

(a) Initial configuration.

(b) Final configuration after 360 cycles.

Figure 17: Density distribution in the femur at baseline and after bone remodeling procedure.

43

(a) Initial densities

(b) K-means clustering initial

(c) Densities after 120 cycles

(d) K-means clustering 120 cycles

(e) Densities after 240 cycles

(f) K-means clustering 240 cycles

Figure 18: Cross section of the femur in the diaphyseal region, 60 mm above the femur’s lower end: (a), (c) and (e) are densities distribution at the initial point, after 120 and 240 cycles, respectively; (b), (d) and (f) correspond to the image segmentation at the initial point, after 120 and 240 cycles, respectively. K-means clustering algorithm was adopted to evaluate the percent increase of the cortical region (defined here as the sum of areas 1 (yellow) and 2 (grey) in (b), (d) and (f)).

44

Table 4: Average ash density (µ), standard deviation (σ) in [g/cm3 ] and percent of each bone type over the procedure’s iterations for the transversal section of the femur illustrated in Fig. 18. Source: Elaborated by the author.

Bone type

Initial configuration

120 cycles 3

240 cycles 3

Percentage

µ ± σ[g/cm ]

Percentage

µ ± σ[g/cm ]

Percentage

µ ± σ[g/cm3 ]

Trabecular

28.11%

0.437 ± 0.098

23.32%

0.502 ± 0.066

17.24%

0.477 ± 0.098

Cortical

71.89%

1.147 ± 0.257

76.67%

1.043 ± 0.179

82.76%

1.245 ± 0.207

Figure 19: Algorithm convergence after approximately 360 iterations. Source: Prepared by the author.

657

or 480 days.

658

Fig. 18 ilustrates a cross section of the femur in the diaphyseal region af-

659

ter convergence. We observe an increase in cortical bone area of approximately

660

7.68% relative to the initial area. The increase was estimated using K-means clus-

661

tering to segment the image and count the number of pixels in each segmented

662

region (adding up areas 1 and 2 shown in Fig. 18). Convergence is shown in

663

Fig. 19 45

(b) Initial densities distribution

(a) Cutting plan

(c) Densities distribution after 240 cycles

Figure 20: Cutting plane in the neck region of the femur and views of densities distribution at the initial and final configurations.

664

Bone remodeling in the neck region of the femur obtained by numerical sim-

665

ulation is illustrated in Fig. 20. Significant change in density distribution and an

666

increase in cortical bone thickness can be observed. The new density distribution

667

depends greatly on the boundary conditions adopted in the simulation.

668

4. DISCUSSION

669

Here we show a mathematical model for bone remodeling that integrates the

670

bone cells populations dynamics and a physiologically based pharmacokinetic

671

model. This is a first attempt to incorporate recent advances in mechanobiology 46

672

concerning the molecular pathways that regulate bone remodeling into a model

673

that can contribute to the development of new antiresorptive or anabolic drugs.

674

In the future, the platform presented here can evolve to precision medicine, al-

675

lowing safer and more effective therapeutic treatments.

676

Here we adopt the Lemaire proposal, which include the RANK-RANKL-OPG

677

pathway, as the core of our model. Lemaire’s model was modified to include

678

the mechanical stimulus evaluated at the microscale and a pharmacoregulator

679

function related to the drug concentration at the bone sites.

680

More sophisticated models of the dynamics of bone cells have been published

681

including other signaling molecules. However, the inclusion of more variables in

682

the model should be pondered with caution since the calibration of parameters

683

can become very difficult or even unfeasible.

684

The changes in bone densities distribution, in the cortical bone area and the

685

increase in the BMD biomarker obtained in the simulations indicate that our in-

686

tegrated model captures reasonably well the expected adaptive behavior of bone

687

during remodeling when subjected to mechanical and pharmacological stimuli.

688

A transversal section at the diaphyseal region after convergence of the algorithm

689

showed an increase in the cortical bone area of 7.68% relative to the initial area.

690

The bone cells populations model includes an anti-resorptive agent, showed a

691

BMD biomarker increase of 3.28%, slightly higher than that obtained in the sim-

692

ulation with mechanical stimuli only. This difference was smaller than expected,

693

suggesting that the implementation is more sensitive to mechanical than phar-

694

macological stimulus. Note that the mechanical stimulus is represented in the

695

model as a mechanoregulator function added to the precursor osteoblasts ex-

696

pression stimulating their differentiation into active osteoblasts. Thus, the study 47

697

indicates that the pharmacological stimuli should be better represented by an

698

increase in the rate of proliferation of precursor osteoblasts, controlled, for in-

699

stance, by the calibration of λ in the mechanoregulator function. Indeed, in a sce-

700

nario of osteoporosis, a decreased sensitivity to mechanical loading and damage

701

is associated to the decline in osteocyte density (Geissler et al., 2015; Frost, 1964).

702

Note that osteocytes and their processes make up a highly sensitive mechanosen-

703

sory system (Cowin et al., 1991) that affects the differentiation of precursors os-

704

teoblasts and consequently the bone formation by active osteoblasts. Future im-

705

plementations could explore other ways to stimulate cell differentiation.

706

An increase in the rate of osteoclasts apoptosis produces a decrease in the rate

707

of differentiation of pre-osteoblasts. Therefore, the anti-resorptive action should

708

not be based solely on the increase in the rate of osteoclastic apoptosis as this is

709

not sufficient to represent an increase in BMD obtained with drug administration.

710

This effect was obtained with a modification in RANKL expression including a

711

regulatory function corresponding to the presence of the drug.

712

According to the parameters adopted in the simulation, the pharmacological

713

stimulus (pharmaco-regulatory function) are normalized between [0,1], regard-

714

less of drug concentration. In real situations, there is a therapeutic range related

715

to the mechanism of action of each drug, in which a maximum response with

716

minimal toxicity is obtained for most patients (Rosenbaum, 2016).

717

Further research should focus on the calibration of the model parameters to

718

better reproduce the behavior of commercially available antiresorptive medica-

719

tions for the treatment of osteoporosis such as alendronate, and risedronate.

720

In the proposed model the choice of boundary conditions strongly influences

721

the cortical formation in the process of bone remodeling due to the mechanical 48

722

stimulus via microscopic strain energy density. Because we adopted a mechan-

723

ical loading less complex than the physiological loading, the cortical thickness

724

did not increase in a strictly concentric manner and presented certain variations,

725

as can be observed in Fig. 18. Thus, it is more plausible to evaluate the increase

726

of cortical area by means of image segmentation than to measure its thickness

727

since the latter varies largely depending on the measurement angle adopted. The

728

increase of cortical bone area was measured through image segmentation using

729

K-means clustering, which allowed a practical quantitative assessment of the

730

procedure’s effectiveness.

731

The mechanoregulation function used the microscopic strain energy density

732

instead of the macroscopic one. As Scheiner et al. (2013) points out, there is a sig-

733

nificant dependence of macroscopic stress state on the current vascular porosity

734

fvas . In the case of close-to-zero porosity, the two energy densities (at the macro

735

and microscale) tend to coincide, while for high porosities, large deviations are

736

found, and the value of macroscopic strain energy density (macroSED) underesti-

737

mates the value of microSED. Also, adopting a stimulus at the microscale makes

738

more sense in this model because inhibition or activation of bone remodeling

739

occur at the cellular level.

740

One limitation of the model concerns the representation of the bone as a

741

nonhomogeneous isotropic material. Although bones present an anisotropic be-

742

havior, the model showed good agreement with the bone’s expected adaptive

743

behavior and captured reasonably well both the cells interactions and molecu-

744

lar pathways that regulate the bone remodeling. In the case of subject-specific

745

FE models, the adoption of a nonhomogeneous isotropic distribution of material

746

properties derived from CT data, as performed here is preferable to average val49

747

ues reported in the literature. This choice leads to more accurate results, since the

748

distribution of stresses in bone is highly related to the distribution of mechanical

749

properties in the bone tissue (Taddei et al., 2004).

750

The PBPK model produced different distributions of drug concentration for

751

each bone type, resembling a real situation. Modeling the long term effects of

752

bone seeking agents as bisphosphonates remains a challenge. The complete elim-

753

ination of bisphosphonate from the body, for instance, is reported to take 10 years

754

on average. Note that the time scale for pharmacological models is defined, in

755

general, in hours, while bone remodeling occurs in days, months and years.

756

As previously mentioned, the PBPK model allows the assessment of differ-

757

ent drug administration regimens in order to optimize treatments, ensure ther-

758

apeutic drug dosages, and minimize possible side effects and liver overload. For

759

instance, oral intake requires a much larger amount of drug to achieve the same

760

concentration obtained with a constant intravenous infusion, which might re-

761

sult in more side effects due metabolisation in the liver. The pharmacological

762

model typically deals with time scales much smaller than the remodeling model

763

time scale (months and years); for this reason we use only the local steady state

764

drug concentration in the expression of pharmacological stimulus. However,

765

it should be noted that local steady-state drug concentration depends on fac-

766

tors such as administration regimen (route, dose frequency), drug characteris-

767

tics (bioavailability, half-life, form of action), and individual patient characteris-

768

tics (weight, height, and health condition). Such elements are not directly con-

769

templated by the remodeling model, but they may be partially contemplated in

770

the physiologically-based pharmacokinetic model. As seen in Section 2.1.1, the

771

proposed pharmacological model provides the temporal evolution of drug con50

772

centration in 4 distinct sub-compartments: cortical and trabecular extravascular

773

bone matrix, and cortical and trabecular vascularized bone surface. A uniform

774

distribution of the drug in each compartment is an inherent assumption in com-

775

partmental models. However, given the extreme morphological and physiologi-

776

cal complexity of the bone, the actual spatial distribution of the drug across the

777

various bones that make up the human skeleton is not uniform. Nonetheless, the

778

authors understand that the fact that local drug concentration can be discrimi-

779

nated in 4 different regions of bone tissue represents a significant advance over

780

current remodeling models, which generally lack this information. Future efforts

781

may move toward exploring local aspects of drug delivery.

782

Results of the bone remodeling simulation agree qualitatively with the ex-

783

perimental results that show an increase in bone mass due to physical exercise

784

and pharmacological treatment. The beneficial effects on BMD associated with

785

exercise vary according to exercise modality (impact and aerobic, resistance, bal-

786

ance and proprioception, whole body vibration, aquatic) and intensity, which

787

renders a quantitative comparison challenging. Ma et al. (2013) found significant

788

and positive effects on the BMD of the femur neck in postmenopausal and per-

789

imenopausal women with walking as a single exercise therapy. Another study

790

(Velez et al., 2008) comparing elderly runners with a sedentary control group

791

found a significantly better whole body BMD for the runner group. Also, a

792

high-intensity resistance and impact training program enhances indices of bone

793

strength and functional performance in postmenopausal women with low bone

794

mass, according to (Watson et al., 2018). Regarding pharmacological treatment,

795

Chailurkit et al. (2003) concluded that treatment with alendronate at a regimen of

796

10 mg daily for postmenopausal women with osteoporosis significantly increased 51

797

BMD at all skeletal sites. According to Liberman et al. (1995), the difference in

798

BMD between women receiving this dosage and those receiving placebo was 5.9

799

0.5 percent in the femoral neck and 7.8 0.6 percent in the trochanter. Fig. 20

800

shows an axial section of the femoral neck and the cortical thickening obtained

801

in the remodeling simulation qualitatively agreeing with the experimental results

802

mentioned above.

803

Here we adopted the steady-state value of drug concentration produced by

804

the PBPK model as an input for the bone remodeling algorithm. Notice also

805

that the concentration of the drug in the bone adopted here is intended only

806

to illustrate an application of the model. Further studies are needed in order to

807

compare our results with real clinical data.

808

5. CONCLUSION

809

Here, we describe a mathematical model for bone remodeling integrating dy-

810

namics of bone cells populations and a physiologically based pharmacokinetics

811

model. The results of the simulations show an increase in cortical bone area,

812

considering a transversal section in the diaphyseal region and an increase in the

813

BMD biomarker evaluated by the system of ordinary differential equations. The

814

model seems to be more sensitive to the mechanical stimulus than to the phar-

815

macological one, suggesting that other signaling paths can be explored in future

816

implementations. Here, the pharmacological stimulus was given by the modifica-

817

tion of the RANKL expression and the increase in osteoclast apoptosis. The study

818

suggests that the differentiation of precursor osteoblasts must also be modified

819

to better reproduce the drug effect (for instance, by modifying the mechanoreg-

820

ulator expression). In conclusion, the model is potentially effective in predicting 52

821

the effects of bone-seeking agents used in bone disease treatments, like alen-

822

dronate, and ziledronate. Further research should focus on the calibration of

823

model parameters to better reproduce the behavior of commercially available

824

antiresorptive medications.

825

Acknowledgements

826

The authors would like to thank Prof. Monica Campiteli of Instituto de F´ısica

827

de S˜ao Carlos of the Universidade de S˜ao Paulo for suggestions and revision of

828

the manuscript.

829

Disclosure statement

830

No conflict of interest is declared.

831

Funding

832

None

833

834

Ethical approval Not required

835

REFERENCES

836

Ara´ujo, D. V., de Oliveira, J. H., Bracco, O. L., 2005. Cost of osteoporotic hip

837

fracture in the brazilian private health care system. Arquivos Brasileiros de

838

Endocrinologia & Metabologia 49 (6), 897–901.

53

839

Bahia, M. T., Hecke, M. B., Mercuri, E. G., 2019. Image-based anatomical recon-

840

struction and pharmaco-mediated bone remodeling model applied to a femur

841

with subtrochanteric fracture: A subject-specific finite element study. Medical

842

Engineering & Physics (69), 58–71.

843

Bahia, M. T., Mercuri, E. G., Hecke, M. B., 2016. Fe analyses of a reconstructed

844

femur with mapping of inhomogeneous isotropic material properties from ct

845

images. VII European Congress on Computational Methods in Applied Sci-

846

ences and Engineering, Crete Island, Greece.

847

Barkaoui, A., 2012. Mod´elisation multi´echelle du comportement m´ecano-

848

biologique de l’os humain: de l’ultrastructure au remodelage osseux. Ph.D.

849

thesis, Universit´e d’Orl´eans.

850

Boutroy, S., Bouxsein, M. L., Munoz, F., Delmas, P. D., 2005. In vivo assessment of

851

trabecular bone microarchitecture by high-resolution peripheral quantitative

852

computed tomography. The Journal of Clinical Endocrinology & Metabolism

853

90 (12), 6508–6515.

854

Boutroy, S., Vilayphiou, N., Roux, J.-P., Delmas, P. D., Blain, H., Chapurlat, R. D.,

855

Chavassieux, P., 2011. Comparison of 2d and 3d bone microarchitecture eval-

856

uation at the femoral neck, among postmenopausal women with hip fracture

857

or hip osteoarthritis. Bone 49 (5), 1055–1061.

858

Brown, R. P., Delp, M. D., Lindstedt, S. L., Rhomberg, L. R., Beliles, R. P.,

859

1997. Physiological parameter values for physiologically based pharmacoki-

860

netic models. Toxicology and industrial health 13 (4), 407–484.

861

Burr, D. B., Allen, M. R., 2013. Basic and applied bone biology. Academic Press. 54

862

Carter, D., Fyhrie, D. P., Whalen, R., 1987. Trabecular bone density and loading

863

history: regulation of connective tissue biology by mechanical energy. Journal

864

of Biomechanics 20 (8), 785789–787794.

865

866

Carter, D. R., Hayes, W. C., 1977. The compressive behavior of bone as a twophase porous structure. The Journal of Bone & Joint Surgery 59 (7), 954–962.

867

Chailurkit, L.-o., Jongjaroenprasert, W., Rungbunnapun, S., Ongphiphadhanakul,

868

B., Sae-Tung, S., Rajatanavin, R., 2003. Effect of alendronate on bone mineral

869

density and bone turnover in thai postmenopausal osteoporosis. Journal of

870

bone and mineral metabolism 21 (6), 421–427.

871

Cooper, D. M., Thomas, C. D. L., Clement, J. G., Turinsky, A. L., Sensen, C. W.,

872

Hallgr´ımsson, B., 2007. Age-dependent change in the 3d structure of cortical

873

porosity at the human femoral midshaft. Bone 40 (4), 957–965.

874

Cowin, S. C., Moss-Salentijn, L., Moss, M. L., 1991. Candidates for the

875

mechanosensory system in bone. Journal of biomechanical engineering

876

113 (2), 191–197.

877

Cremers, S., Sparidans, R., den Hartigh, J., Hamdy, N., Vermeij, P., Papapoulos,

878

S., 2002. A pharmacokinetic and pharmacodynamic model for intravenous bis-

879

phosphonate (pamidronate) in osteoporosis. European journal of clinical phar-

880

macology 57 (12), 883–890.

881

882

Doblar´e, M., Garcıa, J., 2002. Anisotropic bone remodelling model based on a continuum damage-repair theory. Journal of Biomechanics 35 (1), 1–17.

883

Eshelby, J. D., 1957. The determination of the elastic field of an ellipsoidal inclu-

884

sion, and related problems. In: Proceedings of the Royal Society of London A: 55

885

Mathematical, Physical and Engineering Sciences. Vol. 241. The Royal Society,

886

pp. 376–396.

887

´ M., Raffaelli, M. P., Bracco, O. L., Takata, E. T., Reis, F. B., Santili, C., Fortes, E.

888

Lazaretti-Castro, M., 2008. High morbid-mortability and reduced level of os-

889

teoporosis diagnosis among elderly people who had hip fractures in s˜ao paulo

890

city. Arquivos Brasileiros de Endocrinologia & Metabologia 52 (7), 1106–1114.

891

Fritsch, A., Hellmich, C., 2007. fiuniversalfimicrostructural patterns in

892

cortical and trabecular, extracellular and extravascular bone materials:

893

Micromechanics-based prediction of anisotropic elasticity. Journal of Theoret-

894

ical Biology 244 (4), 597–620.

895

Frost, H. M., 1964. Dynamics of bone remodeling. Bone biodynamics 315.

896

Ganghoffer, J.-F., 2012. A contribution to the mechanics and thermodynamics of

897

surface growth. application to bone external remodeling. International Journal

898

of Engineering Science 50 (1), 166–191.

899

Garc´ıa-Aznar, J., Rueberg, T., Doblare, M., 2005. A bone remodelling model cou-

900

pling microdamage growth and repair by 3D BMU-activity. Biomechanics and

901

modeling in mechanobiology 4 (2-3), 147–167.

902

Geissler, J. R., Bajaj, D., Fritton, J. C., 2015. American society of biomechanics

903

journal of biomechanics award 2013: cortical bone tissue mechanical quality

904

and biological mechanisms possibly underlying atypical fractures. Journal of

905

Biomechanics 48 (6), 883–894.

906

907

Gieschke, R., Serafin, D., 2013. Development of innovative drugs via modeling with MATLAB. Springer. 56

908

Goda, I., Ganghoffer, J.-F., Maurice, G., 2016. Combined bone internal and ex-

909

ternal remodeling based on Eshelby stress. International Journal of Solids and

910

Structures 94, 138–157.

911

Gowder, S., 2012. Cell Interaction. IntechOpen.

912

URL https://books.google.es/books?id=8hOeDwAAQBAJ

913

Hambli, R., Boughattas, M. H., Daniel, J.-L., Kourta, A., 2016. Prediction of deno-

914

sumab effects on bone remodeling: A combined pharmacokinetics and finite

915

element modeling. Journal of the mechanical behavior of biomedical materials

916

60, 492–504.

917

Helgason, B., Perilli, E., Schileo, E., Taddei, F., Brynj´olfsson, S., Viceconti, M.,

918

2008. Mathematical relationships between bone density and mechanical prop-

919

erties: a literature review. Clinical Biomechanics 23 (2), 135–146.

920

Hofbauer, L., Kuhne, C., Viereck, V., et al., 2004. The opg/rankl/rank system in

921

metabolic bone diseases. Journal of Musculoskeletal and Neuronal Interactions

922

4 (3), 268.

923

ICRP, 1975. Report of the Task Group on Reference Man. ICRP Publication 23.

924

Pergamon Press, Oxford, international Commission on Radiological Protec-

925

tion.

926

927

Jacobs, C. R., 1994. Numerical simulation of bone adaptation to mechanical loading. Ph.D. thesis, Stanford University, Stanford, EUA.

928

Jones, D., Nolte, H., Schol¨ubbers, J., Turner, E., Veltel, D., 1991. Biochemical signal

929

transduction of mechanical strain in osteoblast-like cells. Biomaterials 12 (2),

930

101–110. 57

931

Kaneko, T. S., Bell, J. S., Pejcic, M. R., Tehranzadeh, J., Keyak, J. H., 2004. Mechan-

932

ical properties, density and quantitative ct scan data of trabecular bone with

933

and without metastases. Journal of Biomechanics 37 (4), 523–530.

934

Kaspar, D., Seidl, W., Neidlinger-Wilke, C., Beck, A., Claes, L., Ignatius, A., 2002.

935

Proliferation of human-derived osteoblast-like cells depends on the cycle num-

936

ber and frequency of uniaxial strain. Journal of Biomechanics 35 (7), 873–880.

937

Keller, T. S., 1994. Predicting the compressive mechanical behavior of bone. Jour-

938

nal of Biomechanics 27 (9), 1159–1168.

939

Komarova, S. V., Smith, R. J., Dixon, S. J., Sims, S. M., Wahl, L. M., 2003. Mathe-

940

matical model predicts a critical role for osteoclast autocrine regulation in the

941

control of bone remodeling. Bone 33 (2), 206–215.

942

943

Kroll, M. H., 2000. Parathyroid hormone temporal effects on bone formation and resorption. Bulletin of mathematical biology 62 (1), 163–188.

944

Lemaire, V., Tobin, F. L., Greller, L. D., Cho, C. R., Suva, L. J., 2004. Modeling the

945

interactions between osteoblast and osteoclast activities in bone remodeling.

946

Journal of theoretical biology 229 (3), 293–309.

947

Liberman, U. A., Weiss, S. R., Br¨oll, J., Minne, H. W., Quan, H., Bell, N. H.,

948

Rodriguez-Portales, J., Downs Jr, R. W., Dequeker, J., Favus, M., et al., 1995. Ef-

949

fect of oral alendronate on bone mineral density and the incidence of fractures

950

in postmenopausal osteoporosis. New England Journal of Medicine 333 (22),

951

1437–1444.

952

Lin, J., Duggan, D., Chen, I.-W., Ellsworth, R., 1991. Physiological disposition of 58

953

alendronate, a potent anti-osteolytic bisphosphonate, in laboratory animals.

954

Drug Metabolism and Disposition 19 (5), 926–932.

955

Lin, J. H., Chen, I., Deluna, F. A., Hichens, M., 1992. Renal handling of alen-

956

dronate in rats. an uncharacterized renal transport system. Drug metabolism

957

and disposition 20 (4), 608–613.

958

Louna, Z., Goda, I., Ganghoffer, J.-F., 2019. Homogenized strain gradient remod-

959

eling model for trabecular bone microstructures. Continuum Mechanics and

960

Thermodynamics, 1–29.

961

Ma, D., Wu, L., He, Z., 2013. Effects of walking on the preservation of bone min-

962

eral density in perimenopausal and postmenopausal women: a systematic re-

963

view and meta-analysis. Menopause 20 (11), 1216–1226.

964

Martin, T., 2004. Paracrine regulation of osteoclast formation and activity: mile-

965

stones in discovery. Journal of Musculoskeletal and Neuronal Interactions 4 (3),

966

243.

967

Masarachia, P., Weinreb, M., Balena, R., Rodan, G., 1996. Comparison of the dis-

968

tribution of 3h-alendronate and 3h-etidronate in rat and mouse bones. Bone

969

19 (3), 281–290.

970

McMullin, T. S., 2005. Drugs. In: Reddy, M., Yang, R., Andersen, M. E., Clewell

971

III III, H. J. (Eds.), Physiologically based pharmacokinetic modeling: science

972

and applications. John Wiley & Sons, New Jersey, Ch. 10, pp. 273–296.

973

974

Miller, P. D., 2009. Denosumab: anti-rankl antibody. Current osteoporosis reports 7 (1), 18–22. 59

975

Morgan, E. F., Bayraktar, H. H., Keaveny, T. M., 2003. Trabecular bone modulus–

976

density relationships depend on anatomic site. Journal of Biomechanics 36 (7),

977

897–904.

978

979

980

981

Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta metallurgica 21 (5), 571–574. Nestorov, I., 2003. Whole body pharmacokinetic models. Clinical pharmacokinetics 42 (10), 883–908.

982

Nestorov, I. A., Aarons, L. J., Arundel, P. A., Rowland, M., 1998. Lumping of

983

whole-body physiologically based pharmacokinetic models. Journal of phar-

984

macokinetics and biopharmaceutics 26 (1), 21–46.

985

NIH, 2001. Nih consensus development panel on osteoporosis prevention, diag-

986

nosis, and therapy, march 7-29, 2000: highlights of the conference. South Med

987

J 94 (6), 569–573.

988

O’Flaherty, E. J., 1991a. Physiologically based models for bone-seeking elements:

989

I. rat skeletal and bone growth. Toxicology and applied pharmacology 111 (2),

990

299–312.

991

O’Flaherty, E. J., 1991b. Physiologically based models for bone-seeking elements:

992

Iii. human skeletal and bone growth. Toxicology and applied pharmacology

993

111 (2), 332–341.

994

Padilla, F., Jenson, F., Bousson, V., Peyrin, F., Laugier, P., 2008. Relationships of

995

trabecular bone structure with quantitative ultrasound parameters: In vitro

996

study on human proximal femur using transmission and backscatter measure-

997

ments. Bone 42 (6), 1193–1202. 60

998

Pastrama, M.-I., Scheiner, S., Pivonka, P., Hellmich, C., 2018. A mathemati-

999

cal multiscale model of bone remodeling, accounting for pore space-specific

1000

mechanosensation. Bone 107, 208–221.

1001

Pertinez, H., Chenel, M., Aarons, L., 2013. A physiologically based pharmacoki-

1002

netic model for strontium exposure in rat. Pharmaceutical research 30 (6),

1003

1536–1552.

1004

1005

Pinheiro, M. d. M., 2008. Mortalidade ap´os fratura por osteoporose. Arquivos Brasileiros de Endocrinologia & Metabologia 52 (7), 1071–1072.

1006

Pinheiro, M. d. M., Eis, S. R., 2010. Epidemiology of osteoporotic fractures in

1007

brazil: what we have and what we need. Arquivos Brasileiros de Endocrinolo-

1008

gia & Metabologia 54 (2), 164–170.

1009

Pivonka, P., Zimak, J., Smith, D. W., Gardiner, B. S., Dunstan, C. R., Sims, N. A.,

1010

Martin, T. J., Mundy, G. R., 2008. Model structure and control of bone remod-

1011

eling: a theoretical study. Bone 43 (2), 249–263.

1012

1013

Porras, A. G., Holland, S. D., Gertz, B. J., 1999. Pharmacokinetics of alendronate. Clinical pharmacokinetics 36 (5), 315–328.

1014

Qin, Q.-H., Wang, Y.-N., 2012. A mathematical model of cortical bone remodel-

1015

ing at cellular level under mechanical stimulus. Acta Mechanica Sinica 28 (6),

1016

1678–1692.

1017

Rieger, R., 2011. Mod´elisation m´ecano-biologique par e´ l´ements finis de l’os

1018

trab´eculaire: des activit´es cellulaires au remodelage osseux. Ph.D. thesis, Uni-

1019

versit´e d’Orl´eans. 61

1020

1021

Rosenbaum, S. E., 2016. Basic pharmacokinetics and pharmacodynamics: An integrated textbook and computer simulations. John Wiley & Sons.

1022

Ross, D. S., Mehta, K., Cabal, A., 2017. Mathematical model of bone remodeling

1023

captures the antiresorptive and anabolic actions of various therapies. Bulletin

1024

of mathematical biology 79 (1), 117–142.

1025

1026

R¨uberg, T., 2004. Computer simulation of adaptive bone remodeling. Ph.D. thesis, Technische Universit¨at Braunschweig, Germany.

1027

Scheiner, S., Pivonka, P., Hellmich, C., 2013. Coupling systems biology with mul-

1028

tiscale mechanics, for computer simulations of bone remodeling. Computer

1029

Methods in Applied Mechanics and Engineering 254, 181–196.

1030

Scheiner, S., Pivonka, P., Hellmich, C., 2016. Poromicromechanics reveals that

1031

physiological bone strains induce osteocyte-stimulating lacunar pressure.

1032

Biomechanics and modeling in mechanobiology 15 (1), 9–28.

1033

Scheiner, S., Pivonka, P., Smith, D., Dunstan, C., Hellmich, C., 2014. Mathe-

1034

matical modeling of postmenopausal osteoporosis and its treatment by the

1035

anti-catabolic drug denosumab. International journal for numerical methods

1036

in biomedical engineering 30 (1), 1–27.

1037

Snyder, S. M., Schneider, E., 1991. Estimation of mechanical properties of cortical

1038

bone by computed tomography. Journal of Orthopaedic Research 9 (3), 422–

1039

431.

1040

1041

Stepensky, D., Kleinberg, L., Hoffman, A., 2003. Bone as an effect compartment. Clinical pharmacokinetics 42 (10), 863–881. 62

1042

Szulc, P., Bouxsein, M. L., 2011. Overview of osteoporosis: Epidemiology and

1043

clinical management. Vertebral Fracture Initiative Resource Document PART

1044

I.

1045

Taddei, F., Pancanti, A., Viceconti, M., 2004. An improved method for the auto-

1046

matic mapping of computed tomography numbers onto finite element models.

1047

Medical engineering & physics 26 (1), 61–69.

1048

Trabelsi, N., Milgrom, C., Yosibash, Z., 2014. Patient-specific fe analyses of

1049

metatarsal bones with inhomogeneous isotropic material properties. Journal

1050

of the mechanical behavior of biomedical materials 29, 177–189.

1051

Trabelsi, N., Yosibash, Z., Milgrom, C., 2009. Validation of subject-specific au-

1052

tomated p-fe analysis of the proximal femur. Journal of Biomechanics 42 (3),

1053

234–241.

1054

Velez, N., Zhang, A., Stone, B., Perera, S., Miller, M., Greenspan, S., 2008. The

1055

effect of moderate impact exercise on skeletal integrity in master athletes. Os-

1056

teoporosis International 19 (10), 1457–1464.

1057

Watson, S. L., Weeks, B. K., Weis, L. J., Harding, A. T., Horan, S. A., Beck, B. R.,

1058

2018. High-intensity resistance and impact training improves bone mineral

1059

density and physical function in postmenopausal women with osteopenia and

1060

osteoporosis: the liftmor randomized controlled trial. Journal of Bone and Min-

1061

eral Research 33 (2), 211–220.

1062

Yokley, K., Tran, H. T., Pekari, K., Rappaport, S., Riihimaki, V., Rothman, N.,

1063

Waidyanatha, S., Schlosser, P. M., 2006. Physiologically-based pharmacoki-

63

1064

netic modeling of benzene in humans: a bayesian approach. Risk Analysis

1065

26 (4), 925–943.

1066

Yosibash, Z., Trabelsi, N., Milgrom, C., 2007. Reliable simulations of the human

1067

proximal femur by high-order finite element analysis validated by experimen-

1068

tal observations. Journal of Biomechanics 40 (16), 3688–3699.

1069

1070

Zaoui, A., 2002. Continuum micromechanics: survey. Journal of Engineering Mechanics 128 (8), 808–816.

64

Highlights:     

Cellular micromechanics and pharmacokinetics integrated in a bone remodeling model Drug concentration at bone sites alters the dynamics of bone cell populations Osteoporosis and its treatment simulated in a mechanobiologic model Physiologically based pharmacokinetics model for bone seeking agents Bone tissue modeled by a permeability rate limited tissue compartment model

Ref: JMBBM_2019_1505 Title: A bone remodeling model governed by cellular micromechanics and physiologically based pharmacokinetics Journal: Journal of the Mechanical Behavior of Biomedical Materials

Declaration of interests X The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: