# PUTTING THE "WHY" BACK INTO BONE "ARCHYTECTURE"

EDITED BY: Phil Salmon, Daniel Chappard and Andrew Anthony Pitsillides PUBLISHED IN: Frontiers in Endocrinology

#### *Frontiers Copyright Statement*

*© Copyright 2007-2016 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

*All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-311-0 DOI 10.3389/978-2-88919-311-0

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

### Dedication to quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

## **PUTTING THE "WHY" BACK INTO BONE "ARCHYTECTURE"**

Topic Editors:

**Phil Salmon,** Bruker-microCT, Belgium **Daniel Chappard,** Institut de Biologie en Santé CHU Angers, France **Andrew Anthony Pitsillides,** Royal Veterinary College, UK

Pattern-shifting swarms or "murmurations" of starlings are an example of subtle dynamics of communication leading to emergent pattern formation.

Photo by Free Wallpaper Point.

A large literature exists on trabecular and cortical bone morphology. The engineering performance of bone, implied from its 3d architecture, is often the endpoint of bone biology experiments, being clinically relevant to bone fracture. How and why does bone travel along its complex spatiotemporal trajectory to acquire its architecture? The question "why" can have two meanings. The first, "teleological - why is an architecture advantageous?" – is the domain of substantial biomechanical research to date. The second, "etiological – how did an architecture come about?" – has received far less attention.

This Frontiers Bone Research topic invited contributions addressing this "etiological why" – what mechanisms can coordinate the activity of bone forming and resorbing cells to produce the observed complex and efficient bone architectures? One mechanism is proposed – chaotic nonlinear pattern formation (NPF) which underlies – in a unifying way – natural structures as disparate as trabecular bone, swarms of birds flying or shoaling fish, island formation, fluid

turbulence and others. At the heart of NPF is the fact that simple rules operating between interacting elements multiplied and repeated many times, lead to complex and structured patterns. This paradigm of growth and form leads to a profound link between bone regulation and its architecture: in bone "the architecture is the regulation". The former is the emergent consequence of the latter.

Whatever mechanism does determine bone's developing architecture has to operate at the level of individual sites of formation and resorption and coupling between the two. This has implications as to how we understand the effect on bone of agents such as gene products or drugs. It may be for instance that the "tuning" of coupling between formation and resorption might be as important as the achievement of enhanced bone volume.

The ten articles that were contributed to this topic were just what we hoped for – a snapshot of leading edge bone biology research which addresses the question of how bone gets its shape. We hope that you find these papers thought-provoking, and that they might stimulate new ideas in the research into bone architecture, growth and adaptation, and how to preserve healthy bone from gestation and childhood until old age.

**Citation:** Salmon, P., Chappard, D., Pitsillides, A. A., eds. (2016). Putting the "Why" Back into Bone "Archytecture". Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-311-0

# Table of Contents


Michael Doube


Gabriel L. Galea, Sion Hannuna, Lee B. Meakin, Peter J. Delisser, Lance E. Lanyon and Joanna S. Price

*45 Early trabecular development in human vertebrae: overproduction, constructive regression, and refinement*

Frank Acquaah, Katharine A. Robson Brown, Farah Ahmed, Nathan Jeffery and Richard L. Abel


Yuchin Wu, Samer Adeeb and Michael R. Doschak

*66 Modalities for visualization of cortical bone remodeling: the past, present, and future*

Kimberly D. Harrison and David M. L. Cooper

*75 3D porous architecture of stacks of a-TCP granules compared with that of trabecular bone: a microCT, vector analysis, and compression study* Daniel Chappard, Lisa Terranova, Romain Mallet and Philippe Mercier

## Editorial: Putting the "Why" Back into Bone "archytecture"

### *Phil Salmon\**

*Bruker-microCT, Kontich, Belgium*

Keywords: bone architecture, morphometry, growth and development, mechanotransduction, nonlinear pattern formation, remodelling, coupling, bone biomaterials

**The Editorial on the Research Topic** 

### **Putting the "Why" Back into Bone "Archytecture"**

This topic asked a question – can we do more than merely describe bone architecture, or analyze its mechanical performance? Can we get insights into the pattern forming mechanisms that generate bone's adaptive architecture? First, I should express gratitude on behalf of myself and co-editors Daniel Chappard and Andy Pitsillides, to the authors who contributed articles. Out of busy academic schedules, they made time to write papers, which provide interesting and novel answers to the above question, from a number of different angles. Also to the reviewers who gave of their time generously, and to the Frontiers staff for their helpful and patient input.

In order to clarify the meaning of this topic with its strange title (containing a deliberate misspelling of "architecture"), I contributed an article (Salmon) "Non-linear pattern formation in bone growth and architecture." More of a stream of consciousness on chaos and pattern than a conventional introduction-methods-results-discussion paper, it threw out some ideas of how chaos-related nonlinear pattern processes may be evident on bone growth and architecture. The paper by Chappard et al. (3D Porous Architecture of Stacks of ß-TCP Granules Compared with That of Trabecular Bone: A microCT, Vector Analysis, and Compression Study) looked at the 3D architecture of an osteogenic scaffold. In terms of the question of how bone remodeling cells behave spatiotemporally within the bone marrow space, this scenario of the bone scaffold could hardly be more appropriate.

The article by Harrison and Cooper (Modalities for Visualization of Cortical Bone Remodeling: The Past, Present, and Future) directly took up the challenge of moving from the "how" to the "why" in the 3D study of bone remodeling units. They provided a review of studies of osteonal remodeling in cortical bone, covering imaging methodologies, both *ex vivo* and *in vivo*, which can image these structures and the osteonal canal networks. The study by Wu et al. (Using Micro-CT Derived Bone Microarchitecture to Analyze Bone Stiffness – A Case Study on Osteoporosis Rat Bone) studied the link between trabecular bone's 3D architecture (as measured by microCT) and its mechanical performance in mechanical tests or finite element analysis (FEA). Erben in his paper "Hypothesis: coupling between resorption and formation in cancellous bone remodeling is a mechanically controlled event," addressed directly the difficult question of how bone remodels on a mechanistic level to achieve a target architecture, which is mechanically adapted to the loads and load directions it experiences.

The study by Acquaah et al. (Early trabecular development in human vertebrae: overproduction, constructive regression, and refinement) was based on very rare and valuable datasets of prenatal, neonatal, and infant human vertebral bones, from a 19th century anatomical collection. This allowed a unique study of the 4D trajectory of growth and bone architecture at this site from 3 months before birth to 2.5 years age, which could prove extremely useful in examining the interaction of genetically determined "baseline" trabecular architecture with an increasing mechanically adaptive

*Edited and Reviewed by: Jonathan H. Tobias,* 

*University of Bristol, UK*

*\*Correspondence: Phil Salmon phil.salmon@bruker.com*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

*Received: 25 January 2016 Accepted: 28 January 2016 Published: 11 February 2016*

#### *Citation:*

*Salmon P (2016) Editorial: Putting the "Why" Back into Bone "Archytecture". Front. Endocrinol. 7:14. doi: 10.3389/fendo.2016.00014*

component. This connects directly with the paper by Erben on how bone remodeling cells respond to these loads.

Developing further the theme of mechanical loading and skeletal architecture, the paper by Galea et al. "Quantification of alterations in cortical bone geometry using site specificity software in mouse models of aging and the responses to ovariectomy and altered loading," came from the lab of Lance Lanyon who is arguably the founder of the *in vivo* study of skeletal biomechanics and load transduction. Over the last half century, Lanyon and his colleagues have pioneered techniques now used world-wide to measure accurately the strains actually experienced in the bones of animal models under controlled loads. One important insight that has come from the work of this group is that the architectural response of a bone to loading must be studied along the bone as a whole, not restricted to a single site such as the midshaft. Galea and his colleagues found that changes in a long bone, in response to loading, disuse, ovariectomy and aging, are very different at different locations, since the "goal" of the bone modeling response is to change its overall shape.

Non-linear pattern formation in bone originates ultimately from the nature of interactions within and between communities of bone remodeling cells – the osteoclasts and osteoblasts. We were fortunate to have a paper contributed by leading authorities on this subject, Sims and Martin. Their paper asked the critical question: "Coupling signals between the osteoclast and osteoblast: how are messages transmitted between these temporary visitors to the bone surface?" In an illuminating focused short review, a detailed picture is given of the osteoclast-to-osteoblast coupling critical to the emergent morphology of bone. It turns out that there is much more to this coupling than the identification of cytokine and signal-receptor links. Microanatomy of the marrow stroma and structures like the resorption "canopy," which transiently lifts over the BMU site like a protective umbrella, also play a role. The review provides a rich source of possible topics

#### **Conflict of Interest Statement:** Phil Salmon is an employee of Bruker-microCT.

*Copyright © 2016 Salmon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or*  for further research into this coupling. The paper's figure 1 is a wonderfully clear annotated cartoon of the resorption to reversal to formation sequence which will likely find its way into many presentations and thesis introduction chapters where bone remodeling essentials are communicated.

The paper by Gao et al. in the group of Janet Henderson in Montreal addressed the bone pathology scoliosis in which an abnormality of bone growth leads to asymmetric curvature of the spinal column ("Micro CT analysis of spine architecture in a mouse model of scoliosis"). This pathology can only be described with reference to 3D topography – a spinal axis which twists in 3D departing from bilateral symmetry. MicroCT analysis of the spines of a genetic mouse model of scoliosis allowed accurate geometric quantification of this disease, providing an important new tool to shed light on its etiology at a mechanistic level. The paper by Doube "The ellipsoid factor for quantification of rods, plates, and intermediate forms in 3D geometries," develops further the topic of quantitative morphometric measurement of bone architecture. Understanding of the complex architecture of especially trabecular bone requires the use of quantitative parameters to characterize important aspects of its architecture. Doube has proposed a new parameter, the ellipsoid factor (EF), which promises to assess plate-rod trabecular architecture, also providing stereological information critical to adaptive load bearing. This parameter is already challenging the existing paradigm of higher trabecular percent volume always meaning more plate-like structure.

In summary, then, the articles contributed were just what we hoped for a snapshot of leading edge bone biology research, which addresses the question of how bone gets its shape.

### AUTHOR CONTRIBUTIONS

The author confirms being the sole contributor of this work and approved it for publication.

*reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Non-linear pattern formation in bone growth and architecture

### **Phil Salmon\***

Bruker-microCT, Kontich, Belgium

#### **Edited by:**

Daniel Chappard, Centre Hospitalier Universitaire d'Angers, France

#### **Reviewed by:**

Jennifer Tickner, University of Western Australia, Australia Jan Josef Stepan, Charles University in Prague, Czech Republic

#### **\*Correspondence:**

Phil Salmon, Bruker-microCT, Kartuizersweg 3B, 2550 Kontich, Belgium e-mail: phil.salmon@bruker.com

The three-dimensional morphology of bone arises through adaptation to its required engineering performance. Genetically and adaptively bone travels along a complex spatiotemporal trajectory to acquire optimal architecture. On a cellular, micro-anatomical scale, what mechanisms coordinate the activity of osteoblasts and osteoclasts to produce complex and efficient bone architectures? One mechanism is examined here – chaotic non-linear pattern formation (NPF) – which underlies in a unifying way natural structures as disparate as trabecular bone, swarms of birds flying, island formation, fluid turbulence, and others. At the heart of NPF is the fact that simple rules operating between interacting elements, andTuring-like interaction between global and local signals, lead to complex and structured patterns.The study of "group intelligence" exhibited by swarming birds or shoaling fish has led to an embodiment of NPF called "particle swarm optimization" (PSO). This theoretical model could be applicable to the behavior of osteoblasts, osteoclasts, and osteocytes, seeing them operating "socially" in response simultaneously to both global and local signals (endocrine, cytokine, mechanical), resulting in their clustered activity at formation and resorption sites. This represents problem-solving by social intelligence, and could potentially add further realism to in silico computer simulation of bone modeling. What insights has NPF provided to bone biology? One example concerns the genetic disorder juvenile Pagets disease or idiopathic hyperphosphatasia, where the anomalous parallel trabecular architecture characteristic of this pathology is consistent with an NPF paradigm by analogy with known experimental NPF systems. Here, coupling or "feedback" between osteoblasts and osteoclasts is the critical element. This NPF paradigm implies a profound link between bone regulation and its architecture: in bone the architecture is the regulation. The former is the emergent consequence of the latter.

**Keywords: bone, morphogenesis, chaos and non-linear dynamics, architecture, developmental biology**

### **INTRODUCTION: 3D ARCHITECTURE IN BONE RESEARCH**

Scientists from many disciplines get drawn into the study of bone architecture. As well as having direct clinical relevance in relation to mechanical competence of the skeleton and fractures, bone structure exerts a fascination – arising perhaps from its functionality, complexity, and even esthetic qualities. A large literature exists on the topic of trabecular and cortical bone architecture. However, the question of how, in the mechanistic developmental sense, bone has the architecture that it does is seldom asked. Why and how does bone travel along its complex spatiotemporal trajectory to acquire its final form? Bone's architecture has generally been taken as a given – the question that follows is "why is a particular architecture advantageous?" This has been the domain of substantial biomechanical research to date. However, the question"in what way did a particular architecture come about?" and has received far less attention.

In the last decade, micro-computed tomography (micro-CT) has had an energizing effect on research into bone biology, by providing a convenient laboratory based method for imaging bone's three-dimensional (3D) architecture non-destructively (1). micro-CT has become a routine tool in the preclinical testing of pharmaceutical agents effects on bone. The center of gravity of research at the leading edge of micro-CT technology within bone research, and also that of related technologies such as synchrotron CT and micro-MRI, has been occupied by the field of bioengineering, and the search for answers as to the precise source of structural or architectural failure leading to broken bones.

The engineering performance of the bone, implied from its 3D architecture, is thus one measured endpoint of bone biology experiments. In general, biologists study bone regulation (involving physiological, genetic, endocrine, and cytokine factors) and the bio-engineers study the architecture. This separation into camps fails to address the direct and profound link that, in fact, exists between bone metabolic regulation and its architecture: in bone the architecture *is* the regulation. The former is the emergent consequence of the latter. How is this?

### **NON-LINEAR PATTERN FORMATION**

For instance, why is trabecular bone trabecular? This is similar to the question – why do large flocks of birds such as starlings sometimes create moving 3D patterns that are mesmerizing and evocative (**Figure 1**)? These spontaneous evolving architectures are referred to as "murmurations." Who or what "tells" each bird where to fly in order to create these beautiful and highly structured

4D patterns? If this were a military pageant, then the answer would be that it is all consciously drilled and rehearsed, micro-managed down to the tiniest movement. But there is no avian "sergeantmajor" barking order at individual citizen starlings. So whence the form and pattern?

What unites these questions – how bone gets its architecture and how the shifting shapes of swarming starlings are generated – is that both exhibit spontaneous pattern-formation associated with the dynamics of chaos. This phenomenon, in which structured shapes appear as if from nowhere in dynamic systems is referred to as non-linear pattern formation (NPF) (2). At the heart of NPF is the fact that quite simple rules of interaction between interacting elements in a dynamic system, multiplied and repeated many times, can lead to highly complex, non-repeating, and structured patterns (3). How one starling in a swarm responds to its neighbor, how bone formation by an osteoblast is stimulated by resorption by a nearby osteoclast and vice versa (**Figure 2**) – the fine scale dynamics of these interaction, multiplied many times over, yield the "emergent" outcome of the evolving shape both of the swarm of starlings and of trabecular bone.

### **COUPLING AND FEEDBACK BETWEEN OSTEOBLASTS AND OSTEOCLASTS**

In the same way that birds in a swarm respond to their neighbors, there is increasing evidence in bone remodeling of active communication both ways between osteoblasts and osteoclasts, better described as coupling or feedback (5). The first major feedback link to be found in bone remodeling was the RANK–RANKL– osteoprotegerin (OPG) system (6). Since that discovery, more

pathways of coupling between osteoblasts and osteoclasts have been uncovered – several of them appear to be connected with the gp130 co-receptor subunit (7). So active feedback between bone formation and resorption is clearly a reality.

This coupling, alternatively referred to as "feedback," is an important component of NPF, and thus the discovery of an increasing number of routes of feedback between bone remodeling cells increases the possibility that NPF may play a role in pattern formation in bone. Some of the ways in which this can happen are explored below.

#### **ROLE OF OSTEOCYTES IN BONE REMODELING**

In recent years, it is also becoming clear that the osteocyte cells, previously seen as passive passengers in remodeling bone, are themselves active agents of signaling and coordination of the activity of osteoblasts and osteoclasts (8). For instance, advanced "4D" micro-CT analysis of remodeling bone by sequential *in vivo* scanning, combined with image coregistration to identify both forming and resorbing regions, and subsequent spatially resolved biochemical analysis, show, for instance, that RANKL expression is increased specifically within a volume of bone that is shortly to undergo resorption [Ref. (9), personal communication]. Likewise, forming regions were found to have elevated OPG expression. This excellent technical work shows that the important bone signaling molecules operate in a precise spatiotemporal context, and can originate from osteocytes within bone.

### **ALAN TURING AND THE "REACTION–DIFFUSION" MODEL**

Alan Turing provided a key insight into the regulation of the development of complex biological tissues by his reaction–diffusion

**FIGURE 2 |The dynamics of interaction and coupling between osteoblasts and osteoclasts give rise to the complex evolving 4-D pattern of trabecular bone, just as the responses between starlings give rise to the highly patterned swarming murmurations**.

model (10), proposed in 1952, essentially an NPF system. The essence of Turing's model was very simple – two biochemical agents operate in a tissue, one promotes cell growth and operates at short range, the second inhibits growth and operates at long range. Turing demonstrated mathematically that complex patterns emerged from the operation of these simple rules. Examples of Turing patterns are shown in **Figure 3**. Turing's insight has proved foundational in biological morphogenesis. Much work since then has confirmed the operation of developed versions of the Turing type reaction–diffusion systems in biological development and morphogenesis (11). For instance, WNT and DKK signaling were shown to represent a Turing reaction–diffusion system in the setting of spacing between murine hair follicles (12). Of course WNT and DKK are well-known players in bone regulation also. Kondo and Miura also referred to the role of "nodal" and "lefty" in left-right asymmetry and that of TGF-beta/FGF and BMPs in tooth pattern, lung branching, and skeletal limb pattern (the *hox* system), as examples of Turing type reaction–diffusion systems.

About a decade after Turing's paper on the reaction–diffusion model, in the early 1960s, the meteorologist Edward Lorenz, working with one of the earliest computers on the simulation of weather systems, discovered again that a set of simple equations, repeated many thousand times, resulted in complex and non-repeating ("non-periodic") patterns (13). He further found that the evolution of these complex systems converged mysteriously toward one or more regions of his parameter "phase space." The term "strange attractor" was later coined for this phenomenon (14). Lorenz' work is often referred to as foundational in the study of emergent pattern from NPF systems.

Turing NPF patterns are visible at the very outset of bone calcification and development. Yochelis et al. showed that the first calcification events, which represent the initiation both of trabecular bone embryonic appearance and also pathological

**FIGURE 3 |Turing patterns in biological development**. From Kondo and Miura (11). Reprinted with permission from AAAS. The image of the popper fish is courtesy of Massimo Boyer (www.edge-of-reef.com).

instances of atopic calcification such as in atherosclerotic plaques, are also characterized by labyrinthine Turing patterns [Ref. (15), **Figure 4A**] representing a reaction–diffusion or promoter– inhibitor morphogenetic system. The same author had previously demonstrated the generation of bone-like labyrinthine patterns with experimental chemical systems (16), such as the classic "BZ" (Belousov–Zhabotinsky) oscillating thin film reaction. The"transverse front instability" or "Ising front" shown in **Figure 4B** bears resemblance to the process of endocortical trabecularization – the thinning of cortical bone by transition to trabecular architecture at the endosteum, which has become a focus of interest recently in research into age-related osteoporosis at cortical sites such as the distal radius (wrist) (17, 18).

### **TURBULENCE**

Once one begins to develop a "feel" for non-linear pattern systems, it is not hard to find them elsewhere in bone structures. For instance, one of the phenomena associated with chaos in fluid dynamics is turbulence. Here is an easy way to see the onset of turbulence. Turn on a tap very slowly. At first, only separated drops of water will fall from the tap. Then a smooth linear flow of water will descend. As you slowly open the tap further, the stream of water will start to wobble and snake from side to side. A little further yet, and the linear stream will break up decisively into a chaotic tumbling cascade. This event is well-known to engineers as the laminar-turbulent transition, or the "turbulent wake transition" (19). Much effort is made by designers to control this phenomenon to reduce wind resistance of cars and improve the flight of aircraft.

**FIGURE 4 | (A)** The labyrinthine pattern of initial calcification, shown by Yochelis et al. (15) in a culture of mesenchymal cells; and **(B)** the progress of a chemical model of non-linear pattern formation – the Belousov– Zhabotinsky reaction – resembling transition from cortical to trabecular bone at an endosteal surface, by the same author (16). ©IOP Publishing and Deutsche Physikalische Gesellschaft. CC BY-NC-SA. Copyright ©2002 Society for Industrial and Applied Mathematics. Reprinted with permission. All rights reserved.

Images of examples of laminar-turbulent fluid transitions are shown in **Figure 5**. Among them is an image of a longitudinal micro-CT section through a distal femur of a juvenile rat. The appearance of the primary spongiosal bone immediately "downstream" of the growth plate, in thin, mostly parallel structures, followed a little further from the growth plate by the more complex-chaotic forms of the mature secondary spongiosal trabecular bone, bears resemblance to a laminar-turbulent transition. It is more than just analogy to describe the newly formed bone at the growth plate in terms of a fluid"flowing"away from the growth plate through the metaphysis as the bone extends. In a sense, it is really a highly viscous and slow moving fluid. The bone's architecture changes hour by hour, day by day, and a parcel of bone formed at the growth plate could be millimeters away from it a week later.

Turbulence represents the full onset of chaos, and the complex 3D structure of turbulent flow can be described as a chaotic cellular structure. This might well be a fitting description of trabecular bone, as the spatiotemporally separated formation and resorption loci transform the bone's architecture in four dimensions, essentially unconstrained in space. The visual likeness of the long bone metaphyseal trabecular bone to a flowing fluid developing turbulence provides a hint as to the underlying pattern-forming process.

### **NPF GETS SOCIAL – PARTICLE SWARM OPTIMIZATION**

We can return to the subject of the spontaneous patterns of the swarming birds (**Figure 1**). The study of these far-from-random evolving 3D shapes of bird swarms and fish shoals, has led to a theoretical model representing a form of NPF, called "particle swarm optimization" or PSO (20). This is related to a body of theory called artificial life or "A-Life." Birds such as starlings which form large flocks are observed to "flock synchronously, change direction suddenly, scatter and re-group iteratively, and finally perch on a target" (21). PSO increases the success rate in foraging for food targets by a bird swarm, or a shoal of fish, and is referred to as "social intelligence" (21). The algorithm of the PSO model simulates quite simple interactions between "particles" (e.g., birds) that include a stochastic (random) component, an adaptive component and positive feedback, in fact, recognizable elements of NPF.

**primary and secondary spongiosal trabecular bone "downstream" of the growth plate in the distal femur of a young rat (E)**.

Without getting into any of the actual maths, the elements of a PSO model can be outlined (21).

### **SOLUTION SPACE**

The solution space is a volume within which a location has to be found that meets criteria for being "optimum." For instance, the volume of air in which the birds are swarming.

### **PARAMETERS**

There are one or more parameters characterizing each location in the solution space.

#### **PROBLEM**

There is a problem to be solved which requires finding the location, within the solution space where one or more parameters have optimal values.

#### **PARTICLE**

The particle exists at a vectorial location within the solution space, and iteratively explores the solution space. The particle senses and registers the value of the parameters at each location that it visits.

### **SWARM**

The swarm is a large group of particles in the solution space. Each particle in the swarm moves at a certain velocity. The particles, as they move around, search for the optimum solution to the problem by referring to previous experiences.

### **PARTICLE'S BEST EXPERIENCE**

Abbreviated to *pbest*, this means the "best" experience of an individual particle, based on the parameter values.

#### **PARTICLE'S NEIGHBOR'S BEST EXPERIENCES**

Abbreviated to *lbest*, or local best, this is the best experience of any of a particle's neighbor particles, up to, for instance, two particles away.

### **SWARM'S BEST EXPERIENCE**

Abbreviated *gbest*, or global best, this is the best experience had by any particle in the swarm.

### **PARTICLE MOVEMENT**

The movement of each particle is partly random, but is also partly influenced by the history of experiences of both that particle and its near neighbors.

#### **ALGORITHM**

The particle moves its position in response to both its own experience and the "swarm intelligence." This is represented by a mathematical weighting toward locations corresponding to both *pbest* and *gbest*.

This PSO model calls to mind a city center on a Friday night full of people looking for the best night club or party venue, comparing experiences of the best DJ's or the worst bouncers, and homing in on the locales where the best time is to be had.

It turns out that the "social influence" or *lbest* is critical to the performance of the swarm (22). Remove *lbest* and the swarm's effectiveness in solving the problem sharply decreases.

The phenomenon of an "attractor" is a feature of chaotic NPF – a structured subset of a system's phase space to which system elements converge. In the chaos literature it is sometimes referred to as a "strange attractor." In PSO models, particle swarms also converge toward an attractor (22). Note that in a NPF system including a PSO model both the phase space and the attractor can be multi-dimensional. Multidimensionality is relevant in a biological cell and tissue context due to the large number of signals and influences that a cell will experience.

Can individual cells, as well as organisms like birds, display "social intelligence" of the PSO type, involving adaptive behavior and interactions and information exchange with neighbors? Evidently they can – bacteria are a well-studied example exhibiting optimized group foraging strategies following PSO or similar models (23–25).

In bone remodeling, could it be that osteoclasts and osteoblasts act as a particle swarm, showing PSO behavior, as they respond to multiple signals operating at different ranges? Could this explain why, for instance, remodeling operates at discreet, concentrated locations – the bone remodeling unit (BMU) elucidated by Frost (26), rather than spreading out more uniformly or randomly on bone surfaces?

Furthermore, it turns out that the 3D architecture of bone marrow is important in cellular participation in bone remodeling. Looking at often-presented signaling cartoons, it is easy to imagine that these processes take place within a uniform featureless soup, or like atoms within a well-mixed gas. Evidence that this is not so was given by de Barros et al. (27) who showed that "migration, proliferation, and differentiation of hematopoietic stem cells (HSCs) are dependent upon a complex 3D bone marrow microenvironment." Thus, as well as endocrine and cytokine actors, there is also 3D landscape involved in bone cellular metabolism. This issue of microenvironment 3D architecture has become recognized as important in the field of regenerative medicine and the fabrication of osteogenic scaffolds for bone tissue repair, where all factors need to be understood which impact on the success of ingrowth of bone progenitor cells into bone repair scaffolds (27, 28).

One can imagine PSO operating between groups of osteoblasts as they negotiate hemopoietic bone marrow environments "looking" individually or collectively for suitable sites to set up a bone modeling unit. Such a paradigm would require evidence that cell populations such as osteoblasts communicate mutually to each other, as well as to other cell populations; such evidence is given by Ziambaras et al. (29); Grellier et al. (30); and Santos et al. (31).

But if we speculate that bone remodeling cells engage in PSO, how could osteocytes, physically locked in their bone lacunae, show swarming behavior? To envision this, we must remember that spaces exist other than the familiar 3D one. Mathematicians use the term "phase space" to describe a system in which each variable parameter is assigned a dimension of its own. Thus, for instance, the expression by osteocytes of RANK, RANKL, OPG, and other signaling molecules could each be assigned a dimension. The full profile of signal expression by an individual osteocyte could be represented by a location in the resulting multi-dimensional phase space (here,one can have as many dimensions as you want). In such a space, our osteocytes are free to soar and swarm like the starlings

discussed earlier, under influence of neighbors, even while locked spatially within their bony lacunae.

### **APPLICATIONS FOR THE NPF PARADIGM IN BONE BIOLOGY**

What has the paradigm of NPF as a driver of bone architecture provided to bone biology, in terms of insights into bone pathologies? One striking example concerns the rare genetic disorder known both as juvenile Pagets disease (JPD) and as Idiopathic Hyperphosphatasia. This consists essentially of a mutation partially or fully inactivating the gene for OPG. OPG is well known as a key component of the RANK-RANKL system of regulation of bone remodeling, OPG being a competitive inhibitor of the signaling to RANK by RANKL. Predictably, withdrawn of this inhibitor or "damper" of remodeling coupling/feedback results in "runaway" remodeling, much increased remodeling rate and severe ensuing osteoporosis (32).

A curious aspect of the bone pathology of JPD shown by Cundy et al. and others by histology is an abnormal parallel arrangement of trabecular plates in iliac crest biopsies, with trabecule aligned in parallel with the iliac cortical wall. Salmon (33) proposed that this parallel trabecular architecture represented the consequence of the removal of damping from resorption–formation coupling or feedback. A physical–chemical analogy to this is the platinumcatalyzed oscillatory oxidation of carbon monoxide on a platinum surface, in which the normally chaotic spiral patterns of oxidation are transformed into an array of parallel lines by the increasing of feedback in an experimentally controlled system (34). This is shown in **Figure 6**. By comparison, 3D micro-CT images of iliac

crest biopsied trabecular bone are shown in **Figure 7** for a JPD patient before and after bisphosphonate treatment (35). The parallel pattern of trabecule is quite striking in the JPD case. In both the JPD trabecular bone and the CO oxidation on platinum, increasing feedback "kills" the chaotic pattern and imposes parallel regularity. Thus, the role of OPG in bone remodeling, as an inhibitor or "damper" of feedback/coupling, can be seen as essential in preserving the normal chaotic architecture of trabecular bone. Its withdrawal can lead to the pathological parallel trabecular pattern in JPD.

Mechanically the parallel trabecular architecture is disastrously weak, as evidenced by the breaking into two halves of the biopsy sample (**Figure 6A**) prior to embedding. Bisphosphonate therapy only partially restored a more normal chaotic trabecular pattern in the central region of the iliac crest, while the trabecule near the cortices remained aligned in parallel with the cortex. Thus, while bisphosphonate slows down the turnover rate (accelerated by JPD) and thereby achieves some limited mitigation of the abnormal architecture, it cannot fully substitute the feedback-damping effect of OPG.

Aside from the exceptional case of JPD, a perspective of NPF establishes a direct link between trabecular architecture and the coupling of remodeling at the level of individual sites of osteoblast formation and osteoclast resorption.

To put it another way, the absence of a pattern-forming perspective in bone research leads to assumptions regarding the link between bone formation and resorption and bone architecture that are inappropriately simplified. An example of this is the belief that increased bone formation rate should lead always to increased trabecular thickness – and inevitable ensuing confusion in interpreting some trabecular morphometric data where this is not the case. The key point is that trabecule cannot be envisioned as static structures, so that the effects of net formation and resorption can simply be added and subtracted arithmetically from surfaces to predict an outcome in thickness. Instead, they are dynamically folding and reshaping their architecture according to emergent non-linear pattern arising from the dynamics of interaction between the agents of bone addition and subtraction – the osteoblasts and osteoclasts. Therefore, data on mere quantitative change to formation and resorption tells you little about the architectural change.

Osteoporosis drugs are often characterized as "anti-resorption" or "anabolic." However, it is impossible to target formation without affecting resorption, and vice versa, since the two are known to be bi-directionally coupled. Through NPF, the spatiotemporal

pattern of change to osteoblast and osteoclast activity and the change to coupling dynamics translates to a complex and not easily predictable change in 3D architecture.

Thus, the "message" given to bone by the gene product or drug or changed mechanical stimulus is not simply a quantitative signal – "get more bone" or "get less bone," but via the agency of NPF becomes a shape or morphological signal – "go to this shape."

In terms of the practice of trabecular bone morphometry, by histology or micro-CT, for instance, the following implications arise from the NPF perspective:


In terms of drug discovery for bone medicine, there are also implications. Principally, it is not necessary to restrict the search of gene targets or drug candidates to ones that will affect bone formation or resorption only quantitatively. More subtle effects on formation–resorption coupling or on the spatiotemporal pattern of osteoblast and/or osteoclast action could also be looked for, which might prove therapeutic even in the absence of a clear quantitative change in bone formation or resorption. Indeed, some studies have found a weakness in the relationship between fracture prevention efficacy of a bone drug and the corresponding reduction in bone mineral density (BMD), a surrogate for the spatial density of trabecular structures (37), and it has been demonstrated that change to bone remodeling rate by itself, independent of the agency of change to BMD, can lower fracture rate (38, 39). Speculation as to the mechanism for this has referred to change to the number of erosion sites ("stress risers") at the bone surface affecting the bone's strength. An alternative possibility, in the light of NPF and the relationship between bone turnover coupling and bone architecture, is that a change in the rate and other characteristics of bone remodeling caused by a therapeutic agent might change the bone's 3D architecture in a way favorable to its mechanical strength. Note that the highly pathological parallel trabecular architecture found in JPD patients, as discussed above, is associated with sharply accelerated bone remodeling.

### **BONE MORPHOLOGY, MECHANICAL LOADING, AND NPF**

The mechanical aspect of bone morphogenesis should not be overlooked of course – bone's architecture is ultimately determined by the demands of loads and torques it is required to resist, according to Wolff's law (40). The primary role of bone is a mechanical one, to provide strength and rigidity and to resist with a margin of safety the mechanical loads encountered in an animal's activities. It has been shown that aspects of bone architecture can be modeled as self-organization in response to mechanical loading (41).

However, if the interactions between formation and resorption include a mechanical term, this will not stop them from being non-linear and thus yielding complex-chaotic emergent pattern when integrated many times over many locations. This helps us understand the difference between nature's designs and our own. To make structures to meet certain engineering requirements we would use H-beams, rectangular plates, round rods, or square meshes, while nature makes chaotic trabecular bone with regional predominant orientations to reflect load directions. We think linear, nature thinks non-linear. In fact, a mechanical component to the "algorithm" determining formation–resorption interactions is likely to be highly non-linear, because all incremental bone growth and remodeling will change the loading environment. And this is the essence of non-linearity that system parameters themselves are changed by the evolution of the system, as was nicely articulated by Gleick (42) in his (highly readable) book "Chaos": "non-linearity means that the act of playing the game has a way of changing the rules."

### **SUMMARY**

To summarize, the idea proposed here is that in bone"the architecture is the regulation." Through the agency of the laws of chaotic NPF, the dynamics of interaction at the smallest level between units of osteoblast formation and osteoclast resorption, multiplied in space and time, give rise to the elaborate patterns of bone, such as trabecular bone, as what is described in chaos terminology as an "emergent" phenomenon. Such pattern phenomena form an important element in developmental biology, and more widely in many spatiotemporal patterns observed in nature. There is also the possibility that cell populations such as osteoblasts, osteoclasts, and osteocytes might act"socially,"responding to mutual signals in order to coordinate their group behavior more optimally, according to the theory of PSO. This perspective needs to be included in the interpretation of bone morphometric observations of the effects of agents of change in bone such as drugs or altered genes. It also focuses attention on the nature of coupling between bone formation and resorption as a fundamental factor directly influencing the architecture of bone.

How could NPF in bone metabolism be tested? The idea of PSO in bone remodeling, for instance, could be investigated by looking for evidence of communication between neighboring cells of the same type – osteoblasts, osteoclasts, and osteocytes. And "*in silico*," computer simulation could aim to reproduce phenomena such as a laminar-turbulent wake transition "downstream" of the growth plate, and transformation between trabecular and cortical bone at endocortical boundaries (trabecular to cortical during bone fetal development and fracture callus formation, the reverse during osteoporosis).

This article falls short of scientific rigor in the sense of logical and experimental proof of this proposal. Instead, a paradigm is proposed, and bone researchers are invited to look for ways in which it might add to the understanding of phenomena of bone architecture observed clinically or experimentally.

### **REFERENCES**


hematopoietic stem cell migration and proliferation in 3D in vitro model. *PloS one* (2010) **5**(2):e9093. doi:10.1371/journal.pone.0009093


**Conflict of Interest Statement:** Phil Salmon is an employee of Bruker-microCT.

*Received: 15 November 2014; paper pending published: 10 December 2014; accepted: 18 December 2014; published online: 20 January 2015.*

*Citation: Salmon P (2015) Non-linear pattern formation in bone growth and architecture. Front. Endocrinol. 5:239. doi: 10.3389/fendo.2014.00239*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology.*

*Copyright © 2015 Salmon. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## The ellipsoid factor for quantification of rods, plates, and intermediate forms in 3D geometries

### **Michael Doube\***

Department of Comparative Biomedical Sciences, The Royal Veterinary College, London, UK

#### **Edited by:**

Phil Salmon, Bruker-microCT, Belgium

#### **Reviewed by:**

Phil Salmon, Bruker-microCT, Belgium Daniel Chappard, Institut de Biologie en Santé CHU Angers, France

#### **\*Correspondence:**

Michael Doube, Department of Comparative Biomedical Sciences, The Royal Veterinary College, Royal College Street, London NW1 0TU, UK e-mail: mdoube@rvc.ac.uk

The ellipsoid factor (EF) is a method for the local determination of the rod- or plate-like nature of porous or spongy continua. EF at a point within a 3D structure is defined as the difference in axis ratios of the greatest ellipsoid that fits inside the structure and that contains the point of interest, and ranges from −1 for strongly oblate (discus-shaped) ellipsoids, to +1 for strongly prolate (javelin-shaped) ellipsoids. For an ellipsoid with axes a ≤ b ≤ c, EF = a/b − b/c. Here, EF is demonstrated in a Java plugin, "Ellipsoid Factor" for ImageJ, distributed in the BoneJ plugin collection. Ellipsoid Factor utilizes an ellipsoid optimization algorithm, which assumes that maximal ellipsoids are centered on the medial axis, then dilates, rotates, and translates slightly each ellipsoid until it cannot increase in volume any further. EF successfully identifies rods, plates, and intermediate structures within trabecular bone, and summarizes the distribution of geometries with an overall EF mean and SD, EF histogram, and Flinn diagram displaying a/b versus b/c. EF is released to the community for testing, use, and improvement.

**Keywords: maximally inscribed ellipsoid, optimization, segmentation, plate, rod,Tb.EF**

### **INTRODUCTION**

The plate- and rod-like shapes observed within trabecular bone may be mechanically important yet represent a challenge to quantification because they are difficult to define and identify in a meaningful way (1–9). The *de facto* standard is the "structure model index" (SMI) which classifies surfaces based on their change in surface area after an infinitesimal dilation, and works well for perfect spheres and cylinders, and other curves of a purely convex nature (9). However, real bone contains a large proportion of concave curvature, and saddle curvature (concave in one direction and convex in the other), which vary as a function of bone volume fraction (BV/TV). In these cases, SMI does not perform well and the final reported figure contains "contamination" from concave portions of the surface. BoneJ (10) implements SMI and uniquely reports the negative contribution of concave surfaces to the SMI sum.

A number of improvements to SMI have been proposed,including Stauber and Müller's volumetric spatial decomposition and Liu et al.'s individual trabecular separation approaches (1, 4). The emphasis of these algorithms is on splitting the structure into discrete elements, then reporting the properties of each element and summing their contribution to the whole. However, except in the most extreme cases, trabecular bone operates as a continuum (a "cellular solid" (11)) rather than as a set of discrete nodes and struts and so the discretization of bony continua into "plates" and "rods"may be artificial as a general solution. In particular, it seems inappropriate when the topology is complex, containing oblique branching and perforation, and when BV/TV rises so that rods and plates are no longer clearly discernible.

### **DEVELOPMENT**

Here, I present an approach for the local classification of threedimensional continua such as trabecular bone based on fitting maximal inscribed ellipsoids. It is intended to overcome limitations of previous concepts and to provide a general solution that makes few assumptions, and that treats bone as a continuum. It is an intuitive extension of the approach used to measure trabecular thickness (Tb.Th), in which the thickness at a point is the diameter of the largest sphere that contains the point and that fits inside the structure (12). The ellipsoid factor (EF) is similarly defined as the difference in axis ratios of the largest ellipsoid, which contains the point and which fits inside the structure. For an ellipsoid (**Figure 1C**) with three semi-axis lengths ("radii") a, b, and c:

$$a \le b \le c, \; EF = \frac{a}{b} - \frac{b}{c} \tag{1}$$

It can be seen from this definition that both *a*/*b* and *b*/*c* have minimum and maximum values of 0 and 1, and so EF can take values from −1 to +1. EF of −1 means that *a*/*b* →0 and *b*/*c* →1: this occurs in extremely oblate (discus-shaped) ellipsoids, which have one short axis and two long axes (**Figure 1D**). EF of +1 means that *a*/*b* →1 and *b*/*c* →0: this occurs in extremely prolate (javelinshaped) ellipsoids, which have one long axis and two short axes (**Figure 1B**). EF of 0 means that *a*/*b* = *b*/*c*, which is an intermediate state between "discus" and "javelin" – spheres satisfy EF = 0 (a:b:c = 1:1:1) along with ellipsoids with axis ratios 1:2:4, 1:3:9, or more generally, a:*q*a:*q* 2 a. A strongly prolate ellipsoid (javelin) would maximally fit a very rod-like feature, while a strongly oblate ellipsoid (discus) would maximally fit a very plate-like feature (**Figure 1A**). The EF at a particular point in the structure is determined by a competition between the overlapping ellipsoids that contain the point. To maintain consistency with the original definition of Local Thickness, here the ellipsoid competition is won by the ellipsoid of greatest volume.

#### **IMPLEMENTATION**

A proof-of-concept implementation has been developed in Java as an ImageJ plugin, and is available, including source code, as part of BoneJ [v1.4.0 and later (10) 1 ]. All code changes are publicly available from GitHub<sup>2</sup> , with developments pushed to the ellipsoid-factor branch, prior to merging with the master branch for general release. Anyone may make improvements to the code, and developers are strongly encouraged to share their changes with the parent project.

Input data are required to be a three-dimensional binary stack. A topology-preserving medial axis thinning is performed with ImageJ's Skeletonize3D plugin (13, 14) and the resulting skeleton points are stored as seed points. For each seed point, a small spherical ellipsoid is instantiated (see BoneJ's Ellipsoid class). A fixed number (default, *n* = 100) of points are drawn on the ellipsoid's surface and the value in the input image at each surface point's coordinate is tested. The ellipsoid is dilated equally along its three axes until at least one surface point lies on a background pixel. When an ellipsoid surface point lies on a background pixel, it is designated as a"contact point."The ellipsoid is then rotated so that one axis aligns with the mean unit vector pointing to the contact point(s), and contracted to lie completely within the foreground. The other two axes are then dilated equally until the number of contact points is at least 1, and the ellipsoid rotated so that a second axis is aligned with the mean vector of contact points.

An iterative cycle of dilation, rotation/translation, and contraction follows. During each iteration, the ellipsoid's volume is determined and if it is the maximal volume found so far, the ellipsoid is stored and used for subsequent iterations. Iteration continues until no further increase in volume has been achieved after a userset number of iterations. Within each iteration, each axis is dilated individually, and the ellipsoid is rotated by a small random rotation [wiggle()], translated [bump()], or turned [turn()] before being contracted until fully fitting within the foreground. Total translation is limited in magnitude so as to prevent ellipsoids drifting far from their seed point, but allowing them to overcome pixelation artifacts relating to discretized 3D space. Translation direction is the mean of the unit vectors from contact points to the center. The turn direction is the torque of unit vector surface normals pushing on the contact points. The general effect of turn() and bump() is for the object boundary to"push back"upon the growing ellipsoid, so that it grows into unfilled space where it is not contacting the boundary. The wiggle() method is included to allow the ellipsoid to overcome ridges and to find growth pathways that might not be found by a strictly analytical approach.

Coordinates outside the image bounds are considered background. Input pixels on the image boundaries tend to seed ellipsoids, which can grow without limit into the out-of-bounds volume. To prevent uncontrolled growth causing a computational halting problem and to avoid results being overly influenced by out-of-bounds space, ellipsoids are culled if they are >50% outside the image bounds; if their volume exceeds the image stack volume; or if their volume has not stabilized before a user-set, large number of iterations.

Following iterative maximal inscribed ellipsoid fitting, an array of all the optimized ellipsoids is formed and sorted in descending order of volume. For each foreground pixel coordinate of the input image, the ellipsoid array is iterated until the first ellipsoid that contains the pixel is found. Because the array is sorted on volume, the first ellipsoid found to contain the point is guaranteed to be the largest. The array index of this ellipsoid is stored in the location of the foreground pixel, forming a 3D map of maximal ellipsoid array indices. Background pixels are set to a large negative number and unmatched foreground pixels for which there is no containing ellipsoid are given the index −1. Further analysis is performed using the ellipsoid array cross-referenced from the 3D map of array indices.

Ellipsoid fitting is substantially more complicated than sphere fitting, due to the additional degrees of freedom. Whereas a sphere is defined by its center and radius alone, an ellipsoid is defined by its center *p*<sup>c</sup> = (*x*c, *y*c,*z* <sup>c</sup>), a 3 × 3 eigenvalue matrix whose diagonal values λ1, λ2, and λ3, relate to the three semi-axis lengths *r* <sup>a</sup>, *r*b, *r* <sup>c</sup> as:

$$r = \frac{1}{\sqrt{\lambda}}\tag{2}$$

<sup>1</sup>http://bonej.org/ef

<sup>2</sup>http://github.com/mdoube/BoneJ/

and a 3 × 3 eigenvector rotation matrix. BoneJ's Ellipsoid class uses the matrix definition of an ellipsoid, rather than the quadratic equation form:

$$ax^2 + by^2 + cz^2 + 2 \text{ dxy} + 2 \text{ fxz} + 2 \text{ gyz} + 2 \text{ hx} + 2 \text{ jy} + 2 \text{ kz} = 1 \text{ (3)}$$

because the matrix form makes translation and rotation transformations trivial to implement and fast to execute. The matrix form also allows fast determination of whether a point lies inside or on an ellipsoid by satisfying the inequality:

$$((X - X\_0)^T H (X - X\_0) \le 1\tag{4}$$

where *H* is the product of the eigenvalue and eigenvector matrices, *X* is the test point, and *X*<sup>0</sup> is the center. To simplify calculation and reduce the number of ellipsoids to optimize, this implementation enforces, like Local Thickness does with spheres (15), that the maximal ellipsoids are centered on the medial axis.

### **RESULTS AND USE**

The Ellipsoid Factor plugin was tested on a Dell T7600 workstation (Dell Products, Bracknell, Berkshire, UK) with 12 CPU cores using three binary stack images, which are available for download at http://bonej.org/ef (**Table 1**, **Figures 2A,E,I**). Two images are X-ray microtomography scans from a previous study (16) that provide a variety of natural geometry. Specifically, trabecular bone from an emu (*Dromaius novaehollandiae*) was selected because it contains large plates separated by well-defined rods (**Figure 2A**), and trabecular bone from a shrew (*Suncus varilla*) was selected to provide a rod-dominated volume (**Figure 2E**). A further synthetic image was constructed by filling hand-drawn ROIs to produce a 3D volume with a small number of intersecting rods and plates (**Figure 2I**). The absolute pixel spacing of these images is irrelevant to EF because being derived from ratios it is a unitless measure; however, the ratio of pixel spacing to feature size is critically important: EF is likely to function decreasingly well as pixel spacing approaches feature size. Mean feature size is reported in **Table 1** as Tb.Th and varies from 14.1 to 16.1 pixels.

Ellipsoid fitting successfully identifies rod- and plate-like regions within trabecular bone, with over 90% of pixels classified by at least one ellipsoid (**Figure 2**). Ridge features, which are neither rod nor plate in the intuitive sense, and which would challenge trabecular separation techniques, are successfully identified. Segmentation of the structure based on EF, or on the individual *a*/*b* and *b*/*c* ratios, is possible, and reveals rod-like features within plates and plate-like regions within rods. Junctions between trabecular elements are identified by intermediate ellipsoids.

Ellipsoid factor can be used as a global measure, like SMI, summarizing the geometry of the whole by calculating the mean pixel value of the EF image. EF forms a much more natural variable than SMI to take the mean of because it varies between −1 and +1 and, in concert with the axis ratio distribution data, the mean ellipsoid can be calculated and displayed directly, unlike SMI, which varies non-linearly from 0 to 4 and for which the resulting mean geometry is difficult to visualize meaningfully. It must be noted that EF, as a simple summary variable, cannot by itself distinguish between a structure formed of equal amounts of rod-like (*a*/*b* →1; *b*/*c* →0, EF→1) and plate-like (*a*/*b* →0; *b*/*c* →1, EF→ −1) ellipsoids, and a structure composed of intermediate ellipsoids (where *a*/*b* = *b*/*c*), because in both conditions the volume-weighted EF of the overall structure tends toward 0. For this reason, an EF histogram should be constructed and interpreted alongside the Flinn diagram, which displays how much of the volume of the structure is described by ellipsoids of particular axis ratios.

The Flinn diagram is more commonly used to model strain and constant volume deformation in geological structures (17, 18), with axis ratios formed by the eigenvalues of the strain tensor. It is a convenient and intuitive method for displaying and analyzing ellipsoid geometry. Here, the Flinn diagram plots the ratios of ellipsoid semi-axis lengths (themselves measured in real spatial units), and for consistency places prolate ellipsoids to the top left and oblate ellipsoids to the lower right.

### **DISCUSSION**

Ellipsoid fitting is much more challenging than sphere fitting, because a sphere is defined by a center and radius,whereas an ellipsoid is defined by a center, three semi-axes, and a rotation. These additional degrees of freedom mean that searching for optimally fitting maximally inscribed ellipsoids is non-trivial in comparison to searching for maximally inscribed spheres, which itself is a computationally intensive task. Here, complexity is reduced by assuming that the maximal ellipsoids are centered on the medial axis, as maximal inscribed spheres are, and so ellipsoids are seeded only from the medial axis and not from every point in the structure.


**Table 1 | Comparison of ellipsoid factor to SMI, BV/TV, and Tb.Th**.

Results of running BoneJ's prototype ellipsoid factor implementation on a 12-core Dell T7600 workstation on three test images. Note the inconsistent relationship between EF and SMI, and the strong negative component to SMI (SMI−), which is nearly 30% of the positive component (SMI+) in the emu image, despite a relatively low volume fraction (BV/TV). The features in these example images are well sampled with a mean thickness of 14.1–16.1 pixels. Processing time increases exponentially for EF (tEF), compared to Tb.Th (tTb.Th), which indicates that improved optimization strategies are required for future EF implementations. All images are available to download from http://bonej.org/ef.

The algorithm described here is a proof of concept that runs sufficiently well to demonstrate the utility of EF in quantifying the rod- and plate-like nature of trabecular bone geometry. However, it must be noted that in its current form it contains some important limitations. The first is speed: the iterative method is relatively slow, requiring tens to hundreds of milliseconds on contemporary hardware to fit each ellipsoid. This is mitigated to some degree by operating in parallel, so that each ellipsoid is optimized independently in a separate thread of the CPU. Large datasets in particular benefit from an approximately linear speedup by increased number of CPU cores. Further speed improvements might be possible by offloading some of the calculations to the GPU, but this may require that the whole image dataset is stored in graphics RAM.

The simple point-probe method used here,where each ellipsoid is tested based on a fixed number of approximately equally spaced points on its surface, has the side-effect of potentially under- or oversampling features in the geometry. This can lead to the optimization ignoring small features such as intratrabecular osteonal canals due to too widely spaced point-probes, or to wasting CPU and memory access cycles due to unnecessarily dense point-probes sampling the same pixel multiple times. An improvement may be to allow the point-probe number and position to vary as a function of ellipsoid size and geometry, so that the inter-point spacing remains similar during ellipsoid dilation and fitting. It may be noted that at the time of its invention, Hildebrand and Rüegsegger's sphere-fitting implementation took around 15 min to process a 286<sup>3</sup> pixel volume on an advanced (for the day) workstation (12). Improved algorithms and hardware mean that a similar Tb.Th measurement job today requires only a few seconds to complete (**Table 1**). I expect that similar improvements in processing speed await EF, should it be found to be useful to the image analysis community.

This implementation, in its current state, does not guarantee complete filling of the input geometry with ellipsoids, unlike the Local Thickness algorithm implemented by Dougherty and Kunzelmann (15). The degree of filling is reported in the log and can be 90% or greater, particularly if all the skeleton points are used to seed ellipsoids and if few or none of the input pixels lie on the image boundaries. An improvement may be to use a seeding strategy other than the current medial axis approach, such as a 3D distance ridge, so that ellipsoids are distributed more evenly throughout the geometry. However, there is still no guarantee that an evenly seeded spherical set of ellipsoids do not optimize away from small surface details due to the optimization strategy, which keeps only the largest discovered ellipsoid centered on the seed point. These residual unclassified pixels might be dealt with by filling them with spheres and giving them an EF of 0, by attaching them to the nearest ellipsoid, or by adding a strategy to the implementation, which aggressively attempts to classify all input pixels by further ellipsoid seeding and fitting.

Finally, the 3D map of EF values is not as smooth as might be expected, which is esthetically disturbing but is an analytically correct result based on the EF definition in which the largest ellipsoid containing a pixel "wins" that pixel. A more balanced and potentially less biased approach might be to allow the nearest ellipsoid to win, or to calculate a distance- or volume-weighted mean at each pixel.

In conclusion, EF is a useful new method for the measurement and segmentation of complicated porous continua such as trabecular bone. It can give a summary of the rod- and plate-like nature of a 3D structure and can identify the dominant geometry at each point within the structure. I suggest the use of the abbreviation Tb.EF to maintain consistency with the standard bone nomenclature (19), when trabecular bone is the input geometry. The current implementation in BoneJ is a working proof-of-concept which the community is encouraged to test and comment on, and upon which improvements in pixel labeling efficiency and computational optimization will be made.

#### **ACKNOWLEDGMENTS**

I thank the two reviewers for their comments, which have helped to improve this manuscript.

#### **REFERENCES**


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 18 December 2014; paper pending published: 15 January 2015; accepted: 26 January 2015; published online: 16 February 2015.*

*Citation: Doube M (2015) The ellipsoid factor for quantification of rods, plates, and intermediate forms in 3D geometries. Front. Endocrinol. 6:15. doi: 10.3389/fendo.2015.00015*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology.*

*Copyright © 2015 Doube. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Micro CT analysis of spine architecture in a mouse model of scoliosis

#### **Chan Gao1,2, Brian P. Chen1,3, Michael B. Sullivan1,2, Jasmine Hui 1,4, Jean A. Ouellet <sup>5</sup> , Janet E. Henderson1,2,5\* and Neil Saran<sup>5</sup>**

<sup>1</sup> Bone Engineering Labs, Research Institute-McGill University Health Centre, Montreal, QC, Canada

<sup>2</sup> Department of Medicine, McGill University, Montreal, QC, Canada

<sup>3</sup> Department of Kinesiology and Physical Education, McGill University, Montreal, QC, Canada

<sup>4</sup> Biotechnology Program, University of British Columbia, Burnaby, BC, Canada

<sup>5</sup> Department of Surgery, McGill University, Montreal, QC, Canada

#### **Edited by:**

Phil Salmon, Bruker-microCT, Belgium

#### **Reviewed by:**

Rob Jurgen van 't Hof, University of Edinburgh, UK Isabel Orriss, Royal Veterinary College, UK

#### **\*Correspondence:**

Janet E. Henderson, Bone Engineering Labs, Montreal General Hospital, McGill University, C9.133, 1650 Cedar Avenue, Montreal, QC H3G 1A4, Canada e-mail: janet.henderson@mcgill.ca

**Objective:** Mice homozygous for targeted deletion of the gene encoding fibroblast growth factor receptor 3 (FGFR3−/−) develop kyphoscoliosis by 2 months of age. The first objective of this study was to use high resolution X-ray to characterize curve progression in vivo and micro CT to quantify spine architecture ex vivo in FGFR3−/<sup>−</sup> mice. The second objective was to determine if slow release of the bone anabolic peptide parathyroid hormone related protein (PTHrP-1-34) from a pellet placed adjacent to the thoracic spine could inhibit progressive kyphoscoliosis.

**Materials and methods:** Pellets loaded with placebo or PTHrP-1-34 were implanted adjacent to the thoracic spine of 1-month-old FGFR3−/<sup>−</sup> mice obtained from in house breeding. X rays were captured at monthly intervals up to 4 months to quantify curve progression using the Cobb method. High resolution post-mortem scans of FGFR3−/<sup>−</sup> and FGFR3+/<sup>+</sup> spines, from C5/6 to L4/5, were captured to evaluate the 3D structure, rotation, and micro-architecture of the affected vertebrae. Un-decalcified and decalcified histology were performed on the apical and adjacent vertebrae of FGFR3−/<sup>−</sup> spines, and the corresponding vertebrae from FGFR3+/<sup>+</sup> spines.

**Results:**The mean Cobb angle was significantly greater at all ages in FGFR3−/<sup>−</sup> mice compared with wild type mice and appeared to stabilize around skeletal maturity at 4 months. 3D reconstructions of the thoracic spine of 4-month-old FGFR3−/<sup>−</sup> mice treated with PTHrP-1-34 revealed correction of left/right asymmetry, vertebral rotation, and lateral displacement compared with mice treated with placebo. Histologic analysis of the apical vertebrae confirmed correction of the asymmetry in PTHrP-1-34 treated mice, in the absence of any change in bone volume, and a significant reduction in the wedging of intervertebral disks (IVD) seen in placebo treated mice.

**Conclusion:** Local treatment of the thoracic spine of juvenile FGFR3−/<sup>−</sup> mice with a bone anabolic agent inhibited progression of scoliosis, but with little impact on kyphosis. The significant improvement in IVD integrity suggests PTHrP-1-34 might also be considered as a therapeutic agent for degenerative disk disorders.

**Keywords: progressive kyphoscoliosis, fibroblast growth factor receptor 3 deficiency, parathyroid hormone-related peptide treatment, radiography, micro-computed tomography, histology**

### **INTRODUCTION**

Scoliosis is a pathological condition characterized by lateral curvature of the spine in the coronal (front–back) plane, whereas kyphosis is excessive angulation of the spine in the sagittal (sideways) plane. Scoliotic disorders are generally classified into congenital forms, such as those arising from fundamental flaws in skeletal patterning in the embryo, those resulting from neuromuscular disorders, and the vast majority of cases that are of undetermined origin (1, 2). The presence of scoliosis before the age of 5 years is defined as early onset scoliosis (EOS) and leads to severe, three dimensional deformation of the spine in ~1/2500 children. The deformity impedes normal lung development and contributes to the primary morbidity associated with the condition (3). While the effects of EOS treatment on morbidity are not fully understood, it has been established that early fusion of the spine to prevent progression of the curve can be further detrimental to pulmonary health (4). Alternative growth-sparing techniques that allow ongoing development of the spine and thoracic cage are therefore encouraged. These techniques rely on bracing, serial casting, and in refractory cases the use of "growing rods" for bio-mechanical correction and inhibition of curve progression. Growing rods are surgical implants placed in the posterior spine to guide longitudinal growth and alignment of the spine. The treatment is fraught with serious complications associated with high financial and psychosocial burden that include significant pain, loss of spinal mobility, and junctional deformity in adjacent vertebrae.

The vast majority of pre-clinical animal models that have been developed to study scoliosis have relied on post-natal interventions including damage to the central nervous system and biomechanical constraints on spine growth (5). However, a genetic basis for some forms of scoliosis is suggested by studies of mice from dams subjected to teratogens like *N*-ethyl-*N*-nitrosourea (6, 7) and carbon monoxide (8), as well as those carrying spontaneous or targeted mutations in genes known to impact the development of neuromuscular and skeletal tissues. For example, mice homozygous for targeted disruption of the gene encoding fibroblast growth factor receptor 3 (*FGFR3* OMIM 134934) exhibit EOS in association with overgrowth of the axial and appendicular skeleton (9, 10). Skeletally mature FGFR3−/<sup>−</sup> mice have also been characterized with reduced cortical bone thickness, defective trabecular bone mineralization, premature joint degeneration, and early arthritis (11, 12).

Parathyroid hormone related protein (PTHrP) was first cloned from malignant tumors removed from patients with a paraneoplastic syndrome known as humoral hypercalcemia of malignancy (13, 14). When released by the tumor into the circulation, the protein, which shares sequence homology with parathyroid hormone (PTH), activates PTH receptors to promote release of calcium from bone and its retention in the kidney (15, 16). Biochemical and physiological analyses determined that PTHrP-1-34 is homologous with PTH-1-34 and binds to the PTH cell surface receptor on target cells, including chondrocytes and osteoblasts (17, 18). A central role for PTHrP in skeletal development was later demonstrated when PTHrP-null mice died at birth with lethal skeletal dysplasia (19) and their heterozygote littermates developed osteopenia as young adults (20). Unlike the catabolic effect mediated by high levels of circulating PTH/PTHrP, the 1–34 peptide mediates anabolic activity on cells when released into the bone micro-environment (21).

FGFR3−/<sup>−</sup> mice have been characterized with osteopenia (11), which has been established as a risk factor for curve progression in patients with scoliosis (22). The objective of this study was to use micro CT to quantify progressive kyphoscoliosis in FGFR3−/<sup>−</sup> mice and to determine if treatment with PTHrP-1-34 could inhibit curve progression. Slow release of the bone anabolic peptide adjacent to the thoracic spine inhibited progression of scoliosis but had little impact on kyphosis or bone volume. The micro CT data was supported by histological analyses of the spine.

### **MATERIALS AND METHODS**

#### **IN VIVO RADIOLOGIC IMAGING OF MOUSE SPINE**

All *in vivo* procedures were conducted as outlined in a protocol approved by McGill Facility Animal Care Committee in compliance with the Canadian Council on Animal Care. Male and female FGFR3+/<sup>+</sup> and FGFR3−/<sup>−</sup> mice used for this study were obtained through in house breeding from a colony maintained for 20 generations on a C3H background. The mice were lightly anesthetized at the indicated times and carefully positioned to capture coronal and sagittal plane radiographs (Kubtek XPERT 80, Milford, CT, USA). Gingko CADx 2.4.1 software (Valladolid, Spain) was used to assess gross anatomical measurements and to measure the magnitude of the largest scoliotic and thoraco-lumbar kyphotic curves, according to the Cobb method (23).

### **TREATMENT WITH SYNTHETIC PTHrP-1-34**

Slow-release pellets measuring 3 mm in diameter and designed to release 1.0 mg PTHrP-1-34 (*N* = 16) or placebo (*N* = 18) over 60 days were manufactured by Innovative Research of America (Sarasota, FL, USA). The pellets were implanted in 4-week-old FGFR3−/<sup>−</sup> mice under general anesthesia on the left side of the spine adjacent to the thoraco-lumbar junction, then maintained for 12 weeks with free access to food and water. Radiographic imaging on anesthetized mice was performed at monthly intervals as described above, and the mice were euthanized by injection of an overdose of anesthetic at 16 weeks of age. Spines were extracted by cutting through the neck and separating the thoracic and lumbar vertebrae from the ribs and pelvis, cleaned of soft tissue, and fixed in 4% paraformaldehyde for 24 h. The fixed spines were then washed in several changes of sterile PBS and stored in PBS at 4°C until the time of micro CT analysis (24).

#### **POST-MORTEM MICRO CT ANALYSES**

High resolution scans of the spine of FGFR3−/<sup>−</sup> (*N* = 16 treated with PTHrP-1-34; *N* = 18 treated with placebo) and FGFR3+/<sup>+</sup> (*N* = 5) mice from C5/6 to L4/5 were captured to evaluate the 3D structure, as well as the rotation of the vertebrae in the curve of FGFR3−/<sup>−</sup> mice. Scans were captured at 8µm spatial resolution on a Skyscan 1172 micro-computed tomograph (micro CT, Bruker, Kontich, Belgium) equipped with a 0.5 mm aluminum filter at a voltage of 55 kV, a current of 180µA, and a power source of 10 W. The rotation step size used for acquiring images was 0.4°. NRecon and CTVol software was used for transverse 2D cross-sectional reconstructions and 3D reconstruction of sub-sets of images.

CTAn software supplied with the instrument was used for quantitative analysis of bone mass and architecture, apical vertebra axial rotation, and geometric properties. Axial rotation of the apical vertebra, as defined by the Scoliosis Research Society (25), was quantified by loading the dataset in Dataviewer and rotating the images to align the anterior and posterior borders of the L5 vertebral body in the horizontal plane. The fifth lumbar vertebra then served as the reference point for measuring the angle of rotation. The properly aligned images were saved as a new dataset and loaded into CTAn to locate the apical vertebra, which is the most rotated vertebra (26). The angle of rotation was measured from the horizontal plane to the straight line from the mid-points of the anterior to posterior surfaces of the apical vertebra.

The L5 vertebral body lays outside of the scoliotic curve in all animals and was therefore used as a reference point to assess bone volume/total volume (BV/TV), which is the primary measure of bone quality. The transverse view image datasets for L5 from FGFR3−/<sup>−</sup> and FGFR3+/<sup>+</sup> mice were used to isolate trabecular bone as the ROI by excluding the vertebral arch and cortical bone. Trabecular bone was then segmented into left, middle, and right portions to calculate BV/TV of the outer (left and right) segments independently. A 3D volume of interest (VOI) that included only the middle 50% of the vertebral body was selected to avoid the endplate region. BV/TV values for bone within the VOI was obtained using CTAn. BV/TV for the apical vertebra in FGFR3−/<sup>−</sup> mice was calculated in the same manner.

The concave/convex geometry of the apical vertebra in FGFR3−/<sup>−</sup> spines was evaluated using cross-sectional images isolated in CTAn and loaded for 3D viewing in Dataviewer. The mid-vertebral body coronal image was saved and the lengths of the concave and convex sides measured from the superior to inferior endplates using CTAn to quantify lateral asymmetry. The same measurements were made on the apical-equivalent vertebra in FGFR3+/<sup>+</sup> mice as a control.

#### **HISTOLOGICAL ANALYSES**

After micro CT analyses, the apical and adjacent vertebrae of FGFR3−/<sup>−</sup> (*N* = 16 treated with PTHrP-1-34; *N* = 18 treated with placebo) spines, and the corresponding vertebrae from the FGFR3+/<sup>+</sup> (*N* = 5) spines, were carefully dissected free from the remaining vertebrae, fixed in 4% paraformaldehyde for 24 h, and embedded at low temperature in poly-methyl-methacrylate (PMMA) as described previously (24). Five micron sections were cut from the polymerized blocks on a Leica microtome and stained with von Kossa/Toluidine blue to distinguish mineralized from soft tissue, with Safranin O to delineate intervertebral disks (IVDs) and vertebral end plates, and for alkaline phosphatase (ALP) and tartrate-resistant acid phosphatase (TRAP) in anabolic and catabolic cells, respectively. Image J software was used to quantify ALP and TRAP activity on the concave vs convex aspects.

#### **STATISTICAL ANALYSES**

A Mann–Whitney *U* test for unpaired samples in the IBM SPSS Statistics (Armonk, NJ, USA) package was used to compare differences in the degrees of scoliosis, kyphosis, axial rotation of apical vertebrae, and concave vs convex height of the apical vertebral body. A one way analysis of variance (ANOVA) with *post hoc* Tukey's test was used for comparison of BV/TV. Differences were considered significant at *p* < 0.05.

### **RESULTS**

Kyphoscoliosis was previously noted in juvenile FGFR3−/<sup>−</sup> mice on a mixed C57Bl6 background but was never characterized (9, 10). The genotype was later transferred onto a C3H background to improve longevity and enable phenotyping of adult mice (11, 12). **Table 1** shows the Cobb angle of thoracic vertebrae measured in large cohorts of FGFR3+/<sup>+</sup> and FGFR3−/<sup>−</sup> mice, aged between 1 and 6 months, that were euthanized in our laboratory. Intragroup variations in spine curvature are remarkably small amongst FGFR3+/<sup>+</sup> mice compared with the large variations in groups of FGFR3−/<sup>−</sup> mice. Beyond 1 month, the mean Cobb angle is significantly greater at all ages in mutant FGFR3−/<sup>−</sup> mice compared with FGFR3+/<sup>+</sup> wild type mice (33.30 ± 19.55° vs 5.06 ± 3.39° at 4 months, *p* < 0.01) and appears to stabilize around 4 months of age when the animals reach skeletal maturity.

#### **PROGRESSION OF SCOLIOSIS IN FGFR3**−**/**<sup>−</sup> **MICE**

To analyze scoliosis, cohorts of FGFR3+/<sup>+</sup> and FGFR3−/<sup>−</sup> mice were anesthetized at monthly intervals from 1 to 4 months to capture high resolution coronal and sagittal X-rays of the spine. Serial radiographs of representative mice, shown in **Figure 1**, reveal an established sigmoid curve at 2 months in the FGFR3−/<sup>−</sup> mouse compared with the FGFR3+/+, which increases in severity to 4 months when the mice were euthanized. As shown in **Figure 2**, continuous release of PTHrP-1-34 from pellets implanted adjacent to the thoracic spine at 1 month of age inhibited progression of scoliosis but did not alter the development of kyphosis in FGFR3−/<sup>−</sup> mice.

#### **MICRO CT ANALYSIS OF THE SCOLIOTIC CURVE IN FGFR3**−**/**<sup>−</sup> **MICE**

In keeping with the clinical assessment of scoliosis, we evaluated the lateral displacement of vertebrae along the horizontal plane and their rotation in the axial plane. Using the L5 vertebra as a reference point outside of the curve, the vertebra with maximal lateral displacement and rotation was identified as the apical



Significantly different from FGFR3<sup>+</sup>/<sup>+</sup> mice \*p < 0.01.

vertebra. The 2D cross-sectional micro CT images of a representative FGFR3−/<sup>−</sup> spine shown in **Figure 3** reveal significant rotation of T8 through T12 in the presence of a relatively small lateral dis-

placement of the spine. 3D reconstructions of the lower thoracic

**to the thoracic spine through a small skin incision**. Coronal (left) and

spine of representative FGFR3−/<sup>−</sup> mice show that treatment with placebo (**Figure 3B**) does not correct the lateral curvature or axial rotation, whereas both are improved in mice treated with pellets that release PTHrP-1-34.

scoliosis but had little impact on kyphosis \*p < 0.03

### **TREATMENT WITH PTHrP-1-34 HAS MINIMAL IMPACT ON VERTEBRAL BONE**

To further characterize the defects in the spine of FGFR3−/<sup>−</sup> mice, the composition and architecture of the apical vertebrae were analysed using quantitative micro CT. The 3D reconstructions shown in the coronal plane (**Figures 4A–C**) in **Figure 4** show the reduction in trabecular bone that characterizes the skeletons of FGFR3−/<sup>−</sup> mice. The coronal images (**Figures 4D– F**) confirm the left/right symmetry in the FGFR3+/<sup>+</sup> vertebra and asymmetry and rotation in the FGFR3−/<sup>−</sup> placebo treated mouse (**Figure 4E**), which are partially corrected in the PTHrP-1-34 treated animal (**Figure 4F**). Quantitative analysis of the L5 vertebrae showed a significant reduction in BV/TV in FGFR3−/<sup>−</sup> mice. Independent analysis of the concave and convex sides of the apical vertebrae showed BV/TV in the concave side to be comparable to that in FGFR3+/<sup>+</sup> vertebrae, while that in the convex side remained the same as the L5 reference. **Table 2** shows significant differences in the quantitative micro CT parameters for the concave and convex sides of the apical vertebrae in FGFR3−/<sup>−</sup> spines compared with no difference in FGFR3+/<sup>+</sup> spines. Betweengroup comparisons confirmed previous observations of an overall increase in the length of vertebrae in FGFR3−/<sup>−</sup> mice, but also a significant reduction of height on the concave side of IVDs, measured indirectly as the distance between the endplates of adjacent vertebrae.

#### **TREATMENT WITH PTHrP-1-34 MAINTAINS IVD INTEGRITY**

Qualitative differences in the vertebral bodies and IVDs in response to placebo and PTHrP-1-34 treatment are shown in **Figure 5**. Low magnification images of Safranin O-stained spines (**Figures 5A,B**) clearly demonstrated the improvement in linearity in PTHrP-1-34 treated mice. Higher magnification images revealed wedging in association with deterioration of the mucoprotein in the nucleolus pulposus and apparent disorganization of the annulus fibrosus in placebo treated (**Figure 5C**) compared with PTHrP-1-34 treated (**Figure 5G**) IVDs. Staining of un-decalcified bones with von Kossa/Toluidine blue (**Figures 5D,H**) showed more bone on the concave than convex side of the apical vertebra in both treatment groups. However, as shown in the Safranin O-stained sections there was a clear improvement in morphology in PTHrp-1-34 compared with placebo treated mice (arrows; **Figure 5D** vs **Figure 5H**). Quantification of ALP (**Figures 5E,I**) and TRAP (**Figures 5F,J**) staining in the vertebral end plates also showed no significant differences between placebo and PTHrP-1- 34 treated animals (ALP 8.0 ± 7.7 vs 17.3 ± 10.9; TRAP 10.9 ± 4.4 vs 15.6 ± 3.5).

### **DISCUSSION**

The goal of this study was to analyze the progression of kyphoscoliosis in osteopenic FGFR3−/<sup>−</sup> mice and to determine if localized treatment with a peptide known to promote bone formation *in vivo* could inhibit progression of the deformity. Micro CT and histological analyses showed that sustained release of PTHrP 1-34 from pellets placed adjacent to the thoracic spine inhibited progression of the scoliotic curve, but had little impact on kyphosis or osteopenia.

A technique for quantification of the lateral displacement of the spine in patients with scoliosis was developed by the orthopedic


Significantly different from FGFR3<sup>−</sup>/<sup>−</sup> ; #p < 0.05; ##,\*\*p < 0.01.

Significantly different from contralateral side; \*p < 0.05; \*\*p < 0.01.

surgeon John Cobb in the late 1940s and remains as the preferred technique for assessment of curve severity (23). Measurement of the Cobb angle on serial antero–posterior radiographs of the spines of growing mice revealed 62 of 63 FGFR3−/<sup>−</sup> mice and no FGFR3+/<sup>+</sup> mice met the criteria for scoliosis at the time of euthanasia between 4 and 6 months of age. Rapid progression of the curve coincided with rapid body growth up to 4 months, when the mice reached skeletal maturity, and slowed thereafter.

The phenotype of FGFR3−/<sup>−</sup> mice overlaps with that of individuals carrying a missense mutation in the tyrosine kinase domain of FGFR3 that results in camptodactyly, tall stature, and hearing loss or CATSHL syndrome (OMIM 610474) (27, 28). Phenotypes that include scoliosis are also seen in humans and mice with mutations in genes encoding proteins involved in the downstream signaling pathways of cell surface receptors like FGFR3. One example is von Recklinghausen Disease or Neurofibromatosis

**FIGURE 5 | PTHrP-1-34 treatment inhibits IVD compression in FGFR3**<sup>−</sup>**/**<sup>−</sup> **mice: After micro CT analyzes representative spines from FGFR3**<sup>−</sup>**/**<sup>−</sup> **mice were de-calcified, embedded in paraffin, and thin sections stained with Safranin O [(A–C,G) orange] to delineate IVD and cartilage endplates of the vertebra or for ALP (E,I) orTRAP (F,J) enzyme activity**. Some spines were embedded un-decalcified in plastic for von Kossa/Toluidine blue staining of bone **(D,H)**. Histologic analysis of the thoracic spine **(A,B)** confirmed the positive influence of PTHrP-1-34 treatment in reducing curvature and axial

rotation shown in the radiology. Higher magnification images of apical vertebrae revealed deterioration and asymmetric compression of the IVD in placebo treated FGFR3<sup>−</sup>/<sup>−</sup> mice [asterix **(C)**], which was not apparent in PTHrP-1-34 treated mice. Von Kossa stained sections revealed a concentration of bone on the concave vs convex side of FGFR3<sup>−</sup>/<sup>−</sup> vertebrae [**(D,H)** black] and confirmed the correction of the IVD phenotype with PTHrP-1-34 treatment [**(D,H)** arrows]. There were no clear differences in the ALP and TRAP activity in adjacent vertebrae [**(E)** vs **(I)** and **(F)** vs **(J)**].

Type 1. This autosomal dominant disorder arises from disruption of the neurofibromin gene (OMIM 613113), which encodes a large cytoplasmic protein that regulates several intracellular signaling molecules (29). Another example is disruption of the protein tyrosine phosphatase non-receptor type 11 gene (*PTPN11* OMIM 176876) encoding the Src-homology 2 phosphatase downstream of FGF receptors (30).

In contrast to the relatively rare incidence of CATSHL syndrome caused by FGFR3 inactivation, activating mutations in FGFR3 that result in achondroplasia (31) and hypochondroplasia (32) occur at an estimated frequency between 1/15,000 and 1/40,000 live births. Affected individuals exhibit kyphoscoliosis in association with short-limbed dwarfism in the growing post-natal skeleton (33). Others less fortunate carry activating mutations that results in a peri-natal lethal form of short-limbed dwarfism called Thanatophoric dysplasia (TD) (34). Like the FGFR3−/<sup>−</sup> mice, these babies show little evidence of kyphoscoliosis at birth. The sporadic mutations that result in TD are distributed throughout the extracellular and intracellular components of FGFR3, which exemplifies the heterogenous nature of the disorder. Targeted mutations in the mouse FGFR3 gene have replicated the phenotype of some human skeletal dysplasias to some extent but have failed to clarify the pathogenesis of scoliotic curve progression (35–37). The consistency with which FGFR3−/<sup>−</sup> mice develop kyphoscoliosis during post-natal growth identify them as a valuable resource to further explore its pathogenesis and potential therapeutic options.

In comparison with placebo, treatment of FGFR3−/<sup>−</sup> mice with PTHrP 1-34 appeared to reduce spine curvature suggesting its activity might be complementary to that of FGF signaling through FGFR3 in skeletally mature mice. Previous work showed compound mutant FGFR3−/−/PTHrP−/<sup>−</sup> mice die in the perinatal period with skeletal abnormalities less severe than those seen in the single mutant PTHrP−/<sup>−</sup> littermates (38). This study concluded that PTHrP regulated the pool of growth plate chondrocytes, while signaling through FGFR3 played a more pronounced role in turnover of cartilage to bone at the chondro–osseous junction. Additional studies in skeletally mature FGFR3−/<sup>−</sup> mice revealed osteopenia arising from defective osteoblast mineralization (11). Taken together this work supports the hypothesis that complementary signaling through PTH1R and FGFR3 pathways regulate cartilage and bone growth and metabolism by targeting different cell populations.

Recombinant PTH 1-34 is the only anabolic medication for the treatment of osteoporosis that has been approved for use in humans (39), although its mechanism of action remains poorly defined. The results of decades of intense investigation has determined that pre-osteoblasts in bone marrow, mature osteoblasts adjacent to the bone surface, and osteocytes embedded in bone matrix all respond to intermittent stimulation with amino-terminal fragments of PTH and PTHrP (17, 40, 41). In the current study, PTHrP-1-34 treatment resulted in blunting of the coronal plane scoliosis in FGFR3−/<sup>−</sup> mice as evidenced by a reduction in the discrepancy between concave and convex lengths of the apical vertebra. Rotation of the apical vertebra was also less pronounced in the mice treated with PTHrP-1-34 compared with those receiving placebo. Although the cause of vertebral rotation remains un-defined, there is general agreement amongst physicians treating patients with scoliosis that the degree of rotation of affected vertebra is increased concomitant with curve severity. These changes resulted in less severe scoliotic deformity and closer approximation to the spinal morphology of WT mice, but with little change in the kyphotic phenotype. One possible explanation is that PTHrP-1-34 released from the pellet in the vicinity of the thoraco-lumbar spine promoted the differentiation of osteogenic cells to strengthen the vertebrae. This conjecture is supported by an apparent increase in ALP activity in the vertebral endplates in our study. Further support comes from a systematic search of the literature that suggested intermittent treatment of osteoporosis with PTH 1-34 could improve both the speed of healing and composition of the tissue in cases of spinal fusion (42).

In contrast to the improved symmetry of the apical vertebra and spine in general in FGFR3−/<sup>−</sup> mice treated with PTHrP-1- 34, there was little evidence of an improvement in BV/TV in the apical vertebra compared with the reference L5 vertebra, which was outside of the curve. The higher BV/TV in the concave compared with convex side of the apical vertebra could have resulted from an anabolic response to increased biomechanical strain on the concave side (43). Alternatively, it could result from collapse of the disorganized trabecular network under increased strain, similar to the spine compression fractures seen in patients with severe osteoporosis (44). Given the known biological action of amino terminal PTH and PTHrP on growth plate and articular chondrocytes, the reduction in curvature could also have resulted from the peptide altering growth and/or metabolism of the cartilaginous IVDs. This hypothesis was supported by the reduction in wedging of the IVDs adjacent to the apical vertebra in FGFR3−/<sup>−</sup> mice treated with PTHrP-1-34 compared with placebo. However, it has been proposed in the literature that IVD wedging occurs secondary to progression of the kyphoscoliotic curve in mice with targeted disruption of the *MECOM* gene, which encodes a transcriptional regulatory protein (45). The lack of apparent impact of PTHrP-1- 34 treatment on BV/TV might also reflect inadequacies in the dose and/or release kinetics of the peptide from the pellet. Suffice to say the precise mechanism by which the positive influence of PTHrP-1-34 in correcting progressive scoliosis in growing FGFR3−/<sup>−</sup> mice awaits additional immunochemical analyses *in vivo* and detailed analysis of intact spines subjected to controlled loading *ex vivo.*

### **AUTHOR CONTRIBUTIONS**

CG: planning, data collection, data analysis, manuscript preparation, and analysis. BC: data collection, data analysis, and manuscript preparation. MS: data collection, data analysis, and manuscript preparation. JH: data collection and data analysis. JO: planning and discussion. JH: planning, data collection, data analysis, manuscript preparation, and analysis. NS: planning, data analysis, manuscript preparation, and analysis.

#### **ACKNOWLEDGMENTS**

The authors thank Dr. Ali Esmaeel for providing clinical relevance during the research design and radiologic imaging, Marco Kniefel for assistance with analysis of radiographs, Huifen Wang for breeding, genotyping, and maintaining the colony of FGFR3 mice,Ailian Li for technical assistance with tissue embedding, sectioning, and staining; and Dominique Behrends for critical reading of the manuscript. This study was made possible with financial support from DePuy/Synthes (NS), the Canadian Institutes of Health Research (JH), the FRSQ-sponsored Reseau de recherché en sante buccodentaire et osseuse (RSBO) (JH), and studentships from RSBO (CG, BC, JH), and MITACS (CG).

### **REFERENCES**


hormone-related peptide (PTHrP) gene. *Genes Dev* (1994) **8**:277–89. doi:10. 1101/gad.8.3.277


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 20 December 2014; paper pending published: 09 February 2015; accepted: 06 March 2015; published online: 19 March 2015.*

*Citation: Gao C, Chen BP, Sullivan MB, Hui J, Ouellet JA, Henderson JE and Saran N (2015) Micro CT analysis of spine architecture in a mouse model of scoliosis. Front. Endocrinol. 6:38. doi: 10.3389/fendo.2015.00038*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology.*

*Copyright © 2015 Gao, Chen, Sullivan, Hui, Ouellet, Henderson and Saran. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

### **Natalie A. Sims\* and T. John Martin**

Department of Medicine, St. Vincent's Institute of Medical Research, St. Vincent's Hospital, The University of Melbourne, Fitzroy, VIC, Australia \*Correspondence: nsims@svi.edu.au

#### **Edited by:**

Phil Salmon, Bruker-microCT, Belgium

#### **Reviewed by:**

Florent Elefteriou, Vanderbilt University, USA Jan Tuckermann, Leibniz Institute for Age Research – Fritz Lipmann Institute, Germany

**Keywords: osteoblast, osteoclast, coupling, osteocyte**

Skeletal mass is regulated by two key activities: bone removal (resorption) by hematopoietic lineage osteoclasts and bone matrix formation by mesenchymal lineage osteoblasts. During adult life, these activities occur sequentially on the same surface: a process termed as remodeling. Tiny packets of bone are removed by osteoclasts and replaced by new bone matrix produced by osteoblasts. This continual renewal process allows repair of mechanical imperfections and calcium homeostasis.

The group of cells responsible for remodeling is termed as the basic multicellular unit (BMU) (1). To maintain bone mass at the same level during adulthood, the bone formed in each BMU must replace precisely the amount removed by resorption within that BMU. This stimulation of osteoblast activity in response to resorption is termed "coupling" (2), and it has long been of interest to understand how these two distinct cell types, on the same bone surface but at different times, could be linked so their activities are equal. The BMU and coupling concepts originally included only osteoclasts and osteoblasts, but over recent years, as more cellular contributors to remodeling have been identified (such as T-cells, macrophages, osteocytes, and precursor populations of osteoblasts and osteoclasts), the number of cells in the BMU has expanded (3–5). So too, more signaling pathways within the BMU have been identified (6, 7). All these signals converge on two cell types: the osteoclast and osteoblast, for only those cells are able of bone resorption and bone formation, respectively.

Osteoclasts and osteoblasts are not present on the bone surface simultaneously; the BMU exists, in different forms, at the same location over approximately 6 months in human bone. Early studies using undecalcified bone histology and timed fluorochrome labeling identified that bone resorption in iliac crest trabecular BMUs of adult human bone takes approximately 3 weeks (8), the formation response 3–4 months (9), and between the two activities there is a poorly understood "reversal phase" (10) of approximately 5 weeks (8). In rodents, the duration of this sequence is compressed, but a time delay between resorption and formation still exists: in rat alveolar bone, the reversal phase lasts for approximately 3.5 days (11). These numbers vary also with site, skeletal health, and treatment (12) and in some conditions, including osteoporosis there is an increased duration, or even arrest, of the reversal phase (13, 14). This review will explore mechanisms by which coupling signals may overcome the time delay between bone resorption and bone formation.

### **THE MAIN CLASSES OF COUPLING FACTORS**

There are four main classes of osteoclastderived signals that may promote bone formation in the BMU: (1) matrix-derived signals released during bone resorption, (2) factors synthesized and secreted by the mature osteoclast, (3) factors expressed on the osteoclast cell membrane, and (4) topographical changes effected by the osteoclast on the bone surface. We have reviewed

these different proposed signals extensively elsewhere (6, 7), but will here focus on how each type of signal might influence bone formation after the reversal phase. As it has been noted in economics, there tends to be a proliferation of putative factors that determine process outcomes until they resemble a "zoo of factors," with only a proportion of them being reliably reproduced in later research (15, 16). Since many coupling factors have only been identified in the past 10 years, most require validation by independent groups of researchers, working in multiple systems. This is particularly true of those factors identified in co-culture studies that disregard the time delay between bone resorption and formation.

### **MATRIX-DERIVED FACTORS**

The bone matrix contains a store of latent growth factors, including transforming growth factor β (TGF-β), bone morphogenetic protein 2 (BMP-2), plateletderived growth factor (PDGF), and the insulin-like growth factors (IGFs) (17– 21). All are deposited by osteoblasts during matrix production then released by osteoclastic activity on the bone surface (20), as well as via plasminogen activators (22, 23) and matrix-metalloproteinases (24). These factors, once released from the matrix, are unlikely to remain within the bone microenvironment for some 5– 8 weeks during the reversal phase until they can influence mature osteoblasts upon their arrival at the bone surface. Rather, their main influences would be to stimulate osteoblast progenitors, including their

recruitment (25–27), migration (26, 28– 30), and differentiation (29, 31, 32).

### **OSTEOCLAST-SECRETED FACTORS**

Osteoclasts also secrete products to promote osteoblast precursor recruitment and differentiation and thereby promote bone formation in the BMU (3). Many such "osteoclast-derived coupling factors" have been identified, and some were validated by studies of genetically altered mice: cardiotrophin-1 (33), sphingosine-1-phosphate, Wnt 10b, BMP-6 (34), CTHRC1 (35), and complement factor 3a (C3a) (36). We have described each of these in detail elsewhere (7); importantly, none are produced exclusively by osteoclasts, and many are also produced by other cell types in the vicinity of the BMU. Both active and inactive osteoclasts produce osteoclast-derived coupling factors (37). This concept is clinically important, since anti-resorptive inhibitors of osteoclast activity rather than osteoclast generation, such as cathepsin K inhibitors, may reduce bone resorption without blocking bone formation, leading potentially to an anabolic effect (38). Considering the time delay between bone resorption and formation, these factors cannot be expected to exist in a stable form at the BMU throughout the reversal phase. These factors stimulate not only the differentiated osteoblast but also, in the case of sphingosine-1 phosphate and cardiotrophin-1, stimulate precursor commitment to the osteoblast lineage (33, 34).

### **OSTEOCLAST MEMBRANE-BOUND FACTORS**

It has also been suggested that osteoclasts interact directly via cell surface regulatory proteins with mature osteoblasts to promote their activity. Factors proposed to act in this cell contact-dependent manner include EphrinB2 (39) and Semaphorin D (40). While plausible *in vitro* such cell contact-dependent mechanisms are problematic at the BMU because osteoclasts and osteoblasts are rarely in direct contact during remodeling. If such mechanisms occur, they may exist between osteoclasts and osteoblast precursors in the bone marrow space, or between osteoclasts and osteoblast-lineage cells in the remodeling canopy (41, 42), an anatomical structure observed above the BMU in human samples (43).

### **WHICH CELL TYPE RESPONDS TO THESE OSTEOCLAST-DERIVED MESSAGES?**

### **OSTEOBLAST PRECURSORS**

Whether matrix-released, osteoclastsecreted, or membrane-bound, the osteoblast precursor seems a more likely target for coupling factors than the mature osteoblast. Indeed, mouse genetic experiments indicate that active TGF-β release during bone resorption induces osteoblast precursor migration to prior resorptive sites (29), thus making them available within the BMU for additional differentiation-promoting signals, such as matrix-derived IGF-1 (27). Osteoblast differentiation would be determined by these subsequent events, which are not well understood and may originate from a wide range of cells near the remodeling surface, including macrophages, T-cells, or the vasculature [reviewed in Ref. (7)]. Indeed, osteoclast-derived coupling factors are released by these other cell populations. For example, TGF-β, Semaphorin 4D, Wnt10b, and BMP-2 are released by T-cells (44–46) and macrophages (47, 48), while sphingosine-1-phosphate and PDGF are released by endothelial cells (49, 50). Nevertheless, the delay between osteoclastor matrix-derived coupling factor release and the arrival of osteoblasts on the bone surface (5–8 weeks) is much longer than the time required for osteoblasts to differentiate *in vitro.* Murine stromal cell lines *in vitro* require only 7 days before alkaline phosphatase production commences (51) and <21 days until they express the genes associated with a matrix-embedded osteocyte (52, 53). Recent lineage tracing experiments suggest the process is even faster *in vivo*, with only 6 days passing before α-SMA precursors are incorporate into the bone matrix as osteocytes (54).

Osteoblast precursors, without the aid of any accessory cell, also sense the size and shape of pits formed by osteoclasts and preferentially form bone matrix in those locations (55). Osteoblast precursors respond to changes in surface topography, whether the change is much larger than the cell itself, as occurs with osteoclastic activity, or very much smaller than the cell (56). Altered nanotopography induces osteoblast filopodia formation followed by cytoskeletal changes involving cell adhesion and differentiation. However, whether

osteoblasts preferentially form matrix on disrupted bone because they truly "sense" a change in physical dimensions of the bone or because they detect a change in surface composition is unclear; both scenarios are plausible. However, it is unlikely that the existence of a resorbed space provides all the information required for osteoblasts to initiate and complete refilling of the space with no input from other cell types. If that were true, remodeling would always be balanced and bone mass constant even in the presence of alterations in intercellular signals.

### **REVERSAL CELLS**

Another cell type that may respond to coupling factors are the reversal cells residing on the bone surface between the bone formation and resorption phases (10). These may provide both a physical and a temporal connection between the osteoclast and the osteoblast. They may respond to osteoclastderived coupling factors (57), and pass the necessary signals to the osteoblast precursors as they move onto the bone surface. Interrogating this possibility remains a challenge, since the reversal cells lack specific identifying features that can be targeted by genetic or pharmacological means.

### **REMODELING CANOPY**

The remodeling canopy is another possible target for coupling factors (57). This anatomical structure consists of osteoblast-lineage cells that lift from the bone surface when osteoclastic resorption initiates the remodeling cycle (43). It has been suggested that the canopy is required for completion of the reversal phase, since biopsies from postmenopausal and glucocorticoid-induced osteoporotic patients exhibit incomplete canopies at sites of reversal phase arrest (58, 59). The canopy could also provide a controlled locale in which osteoblast-lineage cells, osteoclasts, and other contributing marrow cells, may exchange factors and influence precursors provided by the associated vasculature (41, 42). The bone remodeling canopy might serve as a way to keep local coupling factor concentrations sufficiently high to allow mesenchymal stem cell recruitment and to stimulate bone formation, or may limit the cellular contributors to the coupling processes. Defining

the canopy's contribution to the coupling process using genetically altered mouse models is limited because this anatomical structure has not been observed in the mouse, the model used most extensively for defining intercellular signaling pathways involved in bone remodeling.

### **OSTEOCYTES**

Is there a role for the osteocyte in transmitting the messages from osteoclast to osteoblast in the BMU? Osteocytes are osteoblast-lineage cells trapped in the bone matrix during bone formation. They form an extensive interconnected network within a fluid-filled canalicular system (60), which may sense and respond to mechanical strain, as well as paracrine and endocrine signals (61). Osteocytes are present within the BMU throughout remodeling, and may therefore provide a system to transfer information from the resorbing osteoclast to bone surface osteoblasts. While few osteoclast-derived "coupling factors" have been shown to directly influence osteocytes, CT-1 reduces expression of sclerostin (52). CT-1 released by osteoclasts might therefore enter the lacunar–canalicular network, and act on those osteocytes closest to the resorptive site. When mature osteoblasts arrive on the bone surface, if sclerostin expression in the local area is still suppressed, this would allow bone formation to occur in this area. The factors that stimulate osteocytes to promote bone formation at resorbed sites may not be limited to CT-1. The mechanical disturbance caused by resorption itself may also "activate" osteocytes, and stimulate them to provide messages that promote bone formation on the same surface. One key question that remains is whether such a paracrine signaling mechanism could overcome the time delay between resorption and formation. For how long can osteocytes "retain" such information?

A mechanism that would not require a delayed communication system involving osteocytes is their detection and response to changes in microstrain that would occur during the remodeling process. Osteocytes would sense not only the increased strain resulting from weakening of the bone as resorption progresses (62) but would also detect when the strain is relieved as the bone is rebuilt by osteoblasts. Such a strain-based model for coupling was proposed some years ago (63) and now that our understanding of osteocyte signaling has improved, possible mediators are now coming to light. For example, sclerostin may mediate this process: sclerostin is significantly reduced by mechanical loading, and increased during unloading (64). Osteocytes may provide the final refining control to ensure that sufficient bone is formed by osteoblasts, generated in response to messages from osteoclasts either directly or via other cells within the BMU.

### **CONCLUDING COMMENTS**

In conclusion, this suggests three mechanisms by which osteoclast-derived coupling signals may overcome the time delay between bone resorption and formation at the BMU (**Figure 1**). Osteoclast-derived factors (either released from the matrix, secreted from the osteoclast, or expressed on the cell membrane) initiate differentiation of very early osteoblast progenitors, with the level of osteoblast activity and numbers of differentiated cells being refined by other factors released by a range of cell types within the BMU. Second, osteoclast-derived factors may act directly on cells that would transmit further signals to both osteoblast precursors and mature osteoblasts; these transmitting cells could include reversal cells, osteoblast-lineage cells in the remodeling canopy, and osteocytes. Finally, physical changes brought about by the osteoclast would provide osteocytic signals required for initiation and completion of the correct level of matrix production by mature osteoblasts on the bone surface.

#### **REFERENCES**


a kinetic model for matrix and mineral apposition. *Metab Bone Dis Relat Res* (1984) **5**:243–52. doi:10. 1016/0221-8747(84)90065-1


(TGF-beta)-binding protein-1 by osteoclasts. A cellular mechanism for release of TGF-beta from bone matrix. *J Biol Chem* (2002) **277**:21352–60. doi:10.1074/jbc.M111663200


formation. *Biochem Biophys Res Commun* (2008) **366**:483–8. doi:10.1016/j.bbrc.2007.11.168


TGF-beta, PGE2, and PAF. *J Clin Invest* (1998) **101**:890–8. doi:10.1172/JCI1112


compartment canopies, reversal phase arrest, and deficient bone formation in post-menopausal osteoporosis. *Am J Pathol* (2014) **184**:1142–51. doi:10.1016/j.ajpath.2013.12.005


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 16 February 2015; paper pending published: 04 March 2015; accepted: 10 March 2015; published online: 24 March 2015.*

*Citation: Sims NA and Martin TJ (2015) Coupling signals between the osteoclast and osteoblast: how are messages transmitted between these temporary visitors to the bone surface? Front. Endocrinol. 6:41. doi: 10.3389/fendo.2015.00041*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology.*

*Copyright © 2015 Sims and Martin. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Quantification of alterations in cortical bone geometry using site specificity software in mouse models of aging and the responses to ovariectomy and altered loading**

*Gabriel L. Galea <sup>1</sup> \*, Sion Hannuna <sup>2</sup> , Lee B. Meakin <sup>1</sup> , Peter J. Delisser <sup>1</sup> , Lance E. Lanyon <sup>1</sup> and Joanna S. Price <sup>1</sup>*

*<sup>1</sup> School of Veterinary Sciences, University of Bristol, Bristol, UK, <sup>2</sup> Faculty of Engineering, University of Bristol, Bristol, UK*

#### *Edited by:*

*Phil Salmon, Bruker-microCT, Belgium*

#### *Reviewed by:*

*Stefan Alexander Jackowski, University of Saskatchewan, Canada David M. L. Cooper, University of Saskatchewan, Canada*

#### *\*Correspondence:*

*Gabriel L. Galea, School of Veterinary Sciences, University of Bristol, Southwell Street, Bristol BS2 8EJ, UK gabriel.galea@bristol.ac.uk*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

> *Received: 26 January 2015 Paper pending published: 23 February 2015 Accepted: 03 April 2015 Published: 23 April 2015*

#### *Citation:*

*Galea GL, Hannuna S, Meakin LB, Delisser PJ, Lanyon LE and Price JS (2015) Quantification of alterations in cortical bone geometry using site specificity software in mouse models of aging and the responses to ovariectomy and altered loading. Front. Endocrinol. 6:52. doi: 10.3389/fendo.2015.00052* Investigations into the effect of (re)modeling stimuli on cortical bone in rodents normally rely on analysis of changes in bone mass and architecture at a narrow cross-sectional site. However, it is well established that the effects of axial loading produce site-specific changes throughout bones' structure. Non-mechanical influences (e.g., hormones) can be additional to or oppose locally controlled adaptive responses and may have more generalized effects. Tools currently available to study site-specific cortical bone adaptation are limited. Here, we applied novel site specificity software to measure bone mass and architecture at each 1% site along the length of the mouse tibia from standard microcomputed tomography (µCT) images. Resulting measures are directly comparable to those obtained through µCT analysis (*R* <sup>2</sup> *>* 0.96). Site Specificity analysis was used to compare a number of parameters in tibiae from young adult (19-week-old) versus aged (19-month-old) mice; ovariectomized and entire mice; limbs subjected to short periods of axial loading or disuse induced by sciatic neurectomy. Age was associated with uniformly reduced cortical thickness and site-specific decreases in cortical area most apparent in the proximal tibia. Mechanical loading site-specifically increased cortical area and thickness in the proximal tibia. Disuse uniformly decreased cortical thickness and decreased cortical area in the proximal tibia. Ovariectomy uniformly reduced cortical area without altering cortical thickness. Differences in polar moment of inertia between experimental groups were only observed in the proximal tibia. Aging and ovariectomy also altered eccentricity in the distal tibia. In summary, site specificity analysis provides a valuable tool for measuring changes in cortical bone mass and architecture along the entire length of a bone. Changes in the (re)modeling response determined at a single site may not reflect the response at different locations within the same bone.

#### **Keywords: bone, mechanical loading, neurectomy, ovariectomy, aging**

### **Introduction**

Mechanical loading is the primary functional determinant of bone mass and architecture. Resident bone cells' responses to local changes in loading-engendered strains site-specifically fine tune architectural features including curvature (1–3), eccentricity (4), cross-sectional thickness (5, 6), and polar moment of inertia, which is a measurement of bone strength (6, 7). In the young, healthy skeleton, this functional adaptation to loading matches bone's form to its load-bearing function through a homeostatic feedback loop commonly referred to as the "mechanostat" (8, 9). This feedback loop maintains bone strains within a tolerated range during habitual levels of loading. However, the mechanostat appears to fail in humans with age, correlating with reduced availability of estrogen leading to the bone fragility characteristic of osteoporosis (10, 11). Thinning of the load-bearing cortices predisposes to fractures, which in long bones occur preferentially at specific sites including the hip and wrist. It is now well documented that the overall geometry of these high-risk sites, particularly of the femoral neck, determines bone resistance to fracture independently of the absolute mass of bone present (12–14).

Despite regional differences in mass and architecture along their length, long bones must locally fine tune their structure to withstand load-bearing as functional units. At the cellular level, this is achieved through modeling (in which the activity of bone forming osteoblasts is independent of resorbing osteoclasts) and re-modeling (in which osteoblast activity follows osteoclastic resorption). Bone mass and architecture in one region influence the mechanical strain environment in other regions (15), such that structural properties at different sites within the same bone are inter-related. For example, the angle of articulation between the femur and tibia predicts the location of bisphosphonate-associated atypical fractures in the femoral diaphysis in humans (16). Experimental studies in rodents using micro-computed tomography (µCT) and histomorphometry analyses have documented the bone formation response following physiological axial loading and demonstrated that it is site specific. For example, in the ulna loading model, an adaptive response is observed distally but not proximally (17) and axial loading in the mouse tibia elicits a response at 37% of the bone's length from the proximal end, but not distally at 75% of the bone's length (6, 18, 19). Consequently, deductions on the effect of loading based on measurements at a single cortical site cannot be generalized to the rest of the bone. This is particularly relevant given the increasing use of *in vivo* µCT scanning of short sections of bone to provide longitudinal data within the same bone.

Unlike the site-specific targeting of bone (re)modeling following loading, systemic and genetic interventions may have indiscriminate, uniform effects on bone mass and architecture and where load-bearing remains constant may intrinsically alter the strain environment. The long-term (re)modeling outcomes of osteogenic or catabolic interventions will therefore be influenced by their interactions with the mechanostat, either suppressing bone's adaptation to loading such that net bone loss occurs or enhancing the sensitivity/responsiveness of the mechanisms involved in the mechanostat such that a lower level of habitual strain is tolerated (18, 20). Since systemic interventions are additive, synergistic, or oppose the site-specific effects of loading, their outcome is expected to have a degree of regionality. This is clearly demonstrated in our recent report that deletion of *Prkc*α in mice leads to age-, gender-, and site-specific alterations in bone structure under the influence of mechanical load-bearing (21). Documentation of these findings in *Prkc*α *−/−* mice required systematic analyses of bone mass and architecture at numerous sites. However, there are limited tools available for analyzing changes in whole-bone cortical architecture following altered loading or systemic interventions such as ovariectomy.

The primary aim of this study was to generate, validate, apply, and freely disseminate a new site specificity software tool, which allows bias-free and convenient quantification of standard measures of bone mass and architecture at multiple cortical bone sites in the mouse tibia. To demonstrate the efficacy of this software, we produced global maps of changes in bone geometry along the length of the mouse tibia in a number of routinely investigated experimental interventions in bone biology research (aging, axial loading, ovariectomy, and disuse).

### **Materials and Methods**

### **Site Specificity Software for the Analysis of μCT Data**

Mouse tibiae were scanned with high resolution µCT (Bruker, Kontich, Belgium) with a voxel size of 4.8 µm as previously reported (6, 19). Images were reconstructed in NRecon (Bruker, Kontich, Belgium) with the following settings: threshold 1.000–1.160, ring artifacts correction 5, beam hardening correction 25%, smoothing 0. Folders containing all the sequentially labeled cross-sectional images from a single bone are defined as the "in path" for site specificity analysis. A new, empty folder labeled "output" is created inside each "in path" folder. The newly developed site specificity program (Data Sheet 1 in Supplementary Material) is then opened for editing in MatLab (R2012a). The site specificity software "in path" (line 13) is set to the folder containing the cross-sectional µCT images of the bone to be analyzed following the format in the example provided (Table S1 in Supplementary Material, Excel file containing a worked out example and a blank sheet into which raw values generated below can be entered). This program is then run as a standard script in MatLab to identify and analyze the µCT slices corresponding to each 1% of the bone's length.

Measures of area given are in pixels, with each pixel representing the scanning voxel size (in our system each pixel represents 4.8 µm). Site specificity isolates the tibia and calculates the area of bone (cortical area, Ct.Ar) and the enclosed non-bone area (marrow area, Ma.Ar, **Figure 1A**) by segmenting the pixel values using the *k*-means algorithm applied to pixel intensity values where *k* is equal to 2. For pixels representing a mixture of tissue types (partial volume effect), their initial classification depends solely on their similarity to the mean of the two clusters. However, in a subsequent processing stage, pixels disconnected from the main bony component are classified as background using connected component analysis. Summation of Ct.Ar and Ma.Ar provides tissue area within the periosteum (Tt.Ar). In cases in which the cortex is breached such that there is no enclosed space, neither Ma.Ar nor Tt.Ar can be calculated for that individual slice but Ct.Ar is still reliable. This is relatively rare in normal mouse tibiae given vessels and cortical defects rarely run directly through the entire cortex in a single cross-section. Errors increase toward the proximal (trabecular bone) and distal (calcaneus) extremes of the bone such

that only the 81 individual slices 4.8 µm thick at each 1% site between 10 and 90% of the bones length were analyzed.

The site specificity program collates all the measures in an Excel (.csv) file within the "output" folder. In addition, it saves binarized images of the tibia without the fibula in the "output" folder as .gif files, which retain the same pixel size (**Figure 1A**) but have a much smaller file size, facilitating subsequent analysis and storage. These .gif files can be exported to the previously validated and freely available BoneJ program in ImageJ (22, 23) to obtain cross-sectional thickness (Cs.Th), moments of inertia (minimum moment of inertia IMin, maximum moment of inertia Imax, and polar moment of inertia PMI through the summation of IMin and IMax), cross-sectional centroid in the *X* and *Y* direction. BoneJ also calculates Feret's maximum diameter, minimum diameter, and angle by rotating the cross-sectional images to identify the maximum and minimum caliper dimensions. Eccentricity was obtained from the maximum and minimum diameters essentially as calculated by Bruker software (24) based on the length of the major axis (Feret Max from BoneJ) and minor axis (Feret Min from BoneJ).

All bones were oriented such that Feret's angle at the 37% site (arbitrarily chosen) was 0°, i.e., parallel to the *x*-axis such that the *x*-axis is anterior–posterior while the *y*-axis is medial–lateral. Centroid was then used to calculate absolute center of mass deviation from the bone's central axis as a measure of curvature. The central axis was defined as an extrapolated straight line running between 5 and 95% of the bone's length (**Figure 1B**). Calculating the absolute deviation from this line defined within the bone itself accounts for any differences in bone position while being scanned, producing values, which were highly consistent between similar mice in a semi-automated manner but which lose information on the direction of curvature along the axis being investigated.

The outputs of the site specificity program are: Ct.Ar, Ma.Ar, Tt.Ar, and binarized images, which can then be analyzed with BoneJ in ImageJ to provide Cs.Th, PMI, eccentricity, and curvature in the anterior–posterior (*x*) or medial–lateral (*y*) directions. Only cortical bone can be analyzed using site specificity. The site specificity and BoneJ analyses for an entire tibia with the data entry form provided (Table S1 in Supplementary Material) take <10 min on a standard laptop computer.

### **Site Specificity Validation**

In order to ensure that bone structural measurements obtained through site specificity analysis are directly comparable with "conventional" µCT analyses, cortical parameters were calculated at four sites in each of six mice using conventional µCT or site specificity analysis. Conventional µCT analysis was performed using CTan (Bruker, Kontich, Belgium) as previously reported by our group (6, 19). In brief, bones were imaged with an X-ray tube voltage of 49 kV and current 200 µA, with a 0.5 mm aluminum filter. The scanning angular rotation was 180° and the angular increment was 0.6°. The voxel size was 4.8 µm isotropically. Images were reconstructed using a modified Feldkamp algorithm in NRecon and opened for analysis in CTan. The regions of interest calculated based on the length of each bone. A region of 100 µCT slices around each region of interest was selected and the fibula manually excluded (versus automatic exclusion using site specificity). Reconstructed cross-sectional gray scale images were segmented into binary images using adaptive local thresholding with a threshold window of 100–255. The images then underwent sequential modifications to reduce noise and improve continuity and then image series underwent removal of white speckles (*<*100 pixels area) and black speckles (*<*20 voxels). The outer periosteal border was isolated by a shrink wrap function stretching over holes with a diameter greater than 10 pixels. The resulting standard cortical bone measures averaged over 100 slices were compared with the Ct.Ar, Tt.Ar, and Ma.Ar values automatically generated by the site specificity software at the single cross-sectional slice at the desired percentage site.

To test the reproducibility of site specificity analysis, a bone was analyzed five times with complete system restarts between analyses. In addition, Ct.Ar calculated by site specificity was compared to the same parameter calculated by BoneJ.

### **Animal Studies**

To investigate the effect of aging, young adult 18-week-old and aged 19-month-old C57B*\*16 mice were obtained from Charles River (*n* = 6 in each age group).

To study the effect of loading and disuse, C57B/16 mice were also obtained from Charles River. For loading experiments, 17 week-old female mice (*n* = 15) were subjected to axial tibial loading of the right limb three times a week for 2 weeks as previously described (6, 19). In brief, the flexed knee and ankle joints are positioned in concave cups; the upper cup containing the knee is attached to an actuator arm of a loading device and the lower cup to a dynamic load cell. The tibia is held in place by a 0.5 N continuous static pre-load. Forty cycles of dynamic load are superimposed with 10 s rest intervals between each cycle. The protocol for one cycle consists of loading to the target peak load, hold for 0.05 s at the peak load, and unloading back to the 0.5 N pre-load at a load of 13.5 N to engender a peak strain of 2,250 µε on the medial surface of the tibia 37% of its length from the proximal end. The left limbs served as normally loaded internal controls. Disuse was induced as previously described (6, 25) through unilateral sciatic neurectomy (SN) (*n* = 6) of the right limb of 17-week-old female mice. Bones were collected 3 weeks following SN. The left limbs served as normally loaded internal controls.

To investigate the effect of estrogen withdrawal, ovariectomy (OVX) was performed as previously described (18, 21) in 8 week-old mice C57Bl/6 mice bred in house. Bones were collected 10 weeks after ovariectomy and compared to non-ovariectomized littermate controls (*n* = 6).

All procedures complied with the UK Animals (Scientific Procedures) Act 1986 undertaken under project license PPL 30-2829 and were reviewed and approved by the University of Bristol ethics committee (Bristol, UK). All mice were allowed free access to water and a maintenance diet containing 0.75% calcium (EURodent Diet 22%; PMI Nutrition International, LLC, Brentwood, MO, USA) in a 12-h light/dark cycle, with room temperature at 21 *±* 2°C. Mice were housed in groups of up to five animals and all cages contained wood shavings, bedding, and a cardboard tube for environmental enrichment.

### **Statistics**

Statistical comparison was by mixed model analysis in SPSS (PASW Statistics, v.18) with bone site as a fixed categorical parameter, the intervention (aging, loading, disuse, OVX, or control) as a fixed effect, and an intervention by site interaction to determine whether the effect of the intervention was site-specific (i.e., significantly different at different sites) or uniform. For loading and disuse analyses, mouse ID was included as a random effect to account for left and right limbs originating from the same mouse. When the effect of the intervention was significant overall, a *post hoc* Bonferroni correction was applied to identify the individual sites at which the effect was significant at *p <* 0.05.

In order to test the ability to detect differences at *p <* 0.05, a simulated 10% global changes in Ct.Ar were statistically analyzed using mixed models. This magnitude of change was selected because OVX caused a *−*10.0 *±* 0.03% change in Ct.Ar averaged across all sites (range *−*4 to *−*20%), described in the Section "Results." To do this, Ct.Ar values from six mice were increased by 10% across all sites (the "Intervention") and analyzed by mixed model analysis with *post hoc* Bonferroni. The overall *p* value associated with the intervention fixed effect and the *p* values comparing each site were then obtained.

Results are presented as the mean *±* SEM.

### **Results**

### **Validation of the Site Specificity Program**

Ct.Ar, Ma.Ar, and Tt.Ar measures obtained from site specificity and conventional µCT analysis were highly correlated (*R* <sup>2</sup> = 0.99 for each parameter, **Figure 2A**), with slight differences between the two likely due to conventional µCT analysis taking the average of 100 CT slices, whereas site specificity analysis analyses individual slices corresponding to each 1% site. Within each site tested, the site specificity value obtained from a single slice will always be different from that generated by conventional µCT analysis if the parameter in question changes rapidly over the 100 lines analyzed. In areas of rapid change, such as at the tibia–fibula junction, the site specificity represents the actual value for the cross-sectional image corresponding to that level rather than an average of slices around the level of interest. Using the statistical approach described, a simulated intervention effect causing 10% change in Ct.Ar could be detected as significant (*p <* 0.05) with *n* = 4 bones. With a sample size of *n* = 6 (the smallest group size used in the studies presented here), a 10% change was detected as significant in 98.8% of sites analyzed.

Parameters obtained from BoneJ were as previously validated (22) and the eccentricity parameter calculated from these was also highly correlated to values obtained from conventional µCT analysis (*R* <sup>2</sup> = 0.96, not shown). Results obtained using this software can therefore be directly compared to data obtained through conventional µCT analysis. Re-analysis of the same bone produced identical results each time (*R* <sup>2</sup> = 1, *t*-test *p* = 1, not shown). Ct.Ar calculated by site specificity and BoneJ were identical (*R* <sup>2</sup> = 1, *t*-test *p* = 1, not shown).

Site specificity analyses illustrate differences in structural parameters along the length of the bone functional unit. The proximal tibia of young mice has a greater cortical area surrounding a larger medullary area, but a narrower cortical thickness as compared with distal regions of the same bone, such that total tissue area is smaller distally than proximally (**Figures 2B–E**). Absolute centroid deviation from the central axis, as a measure of bone curvature, in the anterior–posterior direction is more marked in the proximal to middle region of the bone corresponding to the tibial crest (**Figures 2F,G**). Curvature in the medial–lateral direction is also most marked in the proximal half of the bone, curving back to the midline by 60% of the bone's length from the proximal end (**Figure 2H**). All patterns of curvature and structural parameters were highly consistent between different mice: for example, the tibia–fibula junction occurred within 2% of the bone's length in all young mice analyzed.

### **Effect of Aging**

Aging was associated with smaller cortical area and larger medullary area, but these differences were site-specific such that the greatest differences were observed in the proximal regions of the bone (**Figures 2B,C**). Conversely, bones from aged mice had lower cortical thickness than those from young mice, but this difference was uniform along the length of the bone (**Figure 2D**). Differences in total tissue area between young and aged were not significant (**Figure 2E**).

Aging was associated with significant differences in deviation in centroid from the central axis. Centroid deviation was greater in aged than young adult bones in the anterior–posterior axis, particularly in the proximal half of the bone (**Figure 2G**). Centroid deviation along the medial–lateral axis was greater in aged bones in both the proximal and distal extremes of the region analyzed (**Figure 2H**).

values obtained from six mice at the four indicated sites using site specificity analysis or conventional µCT analysis illustrating direct comparability between these two methods. Inset: radiograph of a tibia and fibula from a young mouse approximately aligned to the corresponding site values in the graphs below. **(B)** Ct.Ar, **(C)** Ma.Ar, **(D)** Cs.Th, and **(E)** Tt.Ar were calculated for bones of young and aged mice. **(F)** Schematic

### *Y*). **(G)** Absolute anterior–posterior and **(H)** medial–lateral centroid deviation in tibiae from young and aged mice. Points represent the mean *±* SEM, *n* = 15 young, and 6 aged. Dots/bars above the graphs indicate sites of significant difference, *p <* 0.05 following Bonferroni correction. Interactions are the site by aging interactions calculated by mixed models; ns not significant, \*\*\**p <* 0.001.

### **Effect of Additional Loading**

We have previously reported that the region of the mouse tibia, which undergoes the greatest osteogenic response following noninvasive axial loading, is localized to 37% of the bone's length from the proximal end, whereas other studies using similar models have analyzed responses at the mid-diaphysis (26–30). Loading increased cortical thickness and cortical area at both sites, whereas total tissue area was only significantly increased more proximally (**Figures 3A–C**). Although the effect of loading on medullary area was significant overall (mixed model *p* = 0.01), and left versus right comparisons were significant (*p <* 0.05) by *t*-test at all sites between 40 and 70% of the bone's length, none of these differences were statistically significant following Bonferroni correction for multiple comparisons (**Figure 3D**). Thus, a study analyzing the effect of loading on medullary area would conclude this as reduced if the 50% site is analyzed, but that it is unchanged if the 37% site is analyzed (Figure S1 in Supplementary Material).

Two weeks of loading did not significantly alter bone centroid deviation in either the anterior–posterior or medial–lateral directions (**Figures 3E,F**).

### **Effect of Disuse**

Bones' response to changes in loading follows a linear continuum between the low strains association with the bone loss and the higher strains associated with bone gain (6). However, unlike loading, SN-induced disuse resulted in uniform loss of cortical thickness throughout the tibial cortex (**Figure 4A**). The reduction

**FIGURE 3 | Analysis of response of the mouse tibia to loading at each site along the bone's length**. Site specificity analysis of **(A)** Cs.Th, **(B)** Ct.Ar, **(C)** Tt.Ar, **(D)** Ma.Ar, **(E)** absolute anterior–posterior, and **(F)** medial–lateral centroid deviation calculated from control left tibiae and

loaded right tibiae from the same mice, *n* = 15. Dots/bars above the graphs indicate sites of significant difference, *p <* 0.05 following Bonferroni correction. Interactions are the site by aging interactions calculated by mixed models; \**p <* 0.05, \*\*\**p <* 0.001.

**FIGURE 4 | Analysis of response of the mouse tibia to neurectomy-induced disuse at each site along the bone's length**. Site specificity analysis of **(A)** Cs.Th, **(B)** Ct.Ar, **(C)** Tt.Ar, **(D)** Ma.Ar, **(E)** absolute anterior–posterior, and **(F)** medial–lateral centroid deviation calculated from

control left tibiae and disused right tibiae from the same mice, *n* = 6. Dots/bars above the graphs indicate sites of significant difference, *p <* 0.05 following Bonferroni correction. Interactions are the site by neurectomy interactions calculated by mixed models; ns not significant, \*\*\**p <* 0.001.

in bone area due to disuse was site specific; greater in the more proximal than distal regions of the bone (**Figure 4B**). Disuse did not alter total tissue area but significantly increased marrow area in a site-specific manner (**Figures 4C,D**). Curvature in the anterior–posterior axis was not significantly altered by disuse whereas medial–lateral curvature was significantly increased at the proximal and distal extremes of the region analyzed, but the site by neurectomy interaction was not statistically significant (**Figures 4E,F**).

### **Effect of Ovariectomy**

Ovariectomy during growth did not alter average cortical thickness achieved in adulthood (**Figure 5A**) but resulted in smaller cortical bone area in a uniform manner (**Figure 5B**). This involved smaller total tissue area and larger medullar cavity in a nonsite-specific manner (**Figures 5C,D**). Ovariectomy did not alter the absolute centroid deviation in the anterior–posterior axis (**Figure 5E**), but significantly increased medial–lateral centroid deviation between approximately 30 and 60% of the bone's length while decreasing it between approximately 65 and 80% of the bone's length from the proximal end (**Figure 5F**). This pattern of change in bone curvature is distinct from that observed following aging or disuse.

### **Site-Specific Influences on Polar Moment of Inertia and Eccentricity**

Aging, axial loading, and ovariectomy all influence overall strength of the bone functional unit. Moments of inertia are structural parameters related to bone strength, which can be calculated from CT images (31). Within each tibia, the polar moment of inertia was greater in the proximal region of the bone than the distal regions beyond the tibia–fibula junction (**Figures 6A–C**). Loading increased polar moment of inertia in a site-specific manner, with significant changes proximally but not distally (**Figures 6A–F**). Tibiae from aged mice followed a similar pattern, with smaller polar moment of inertia in the proximal but not distal regions of the bone relative to young mice (**Figures 6G–I**).

Unexpectedly, ovariectomy also decreased polar moment of inertia in a site-specific manner (**Figures 6J–L**), despite having caused a uniform loss in bone area. We hypothesized that the sitespecific effects on polar moment may be due to uniform bone loss from different baselines along the length of the bone functional unit. Simulated reductions in bone mass at the periosteal or endosteal surface cause relatively greater reductions in polar moment of inertia at the proximal 37% site than at the distal 75% site (Figure S2 in Supplementary Material), indicating that the polar moment of inertia in the distal tibia is relatively unaffected by changes in cortical area compared with the less circular proximal tibia. This suggests that the uniform reduction in bone area caused by ovariectomy would reduce PMI in the proximal tibia to a greater extent than in the distal tibia.

Given moments of inertia influenced by shape, we also investigated whether the interventions tested altered bone shape by calculating the eccentricity parameter, which is a load-bearing responsive geometric parameter (4). Eccentricity was not significantly altered by loading or disuse (not shown), but was significantly affected by both ovariectomy and aging in a sitespecific manner predominantly in the distal regions of the bone (**Figures 7A–C**).

## **Discussion**

Bone mass and architecture are the outcome of bone's functional adaptation to local load-bearing occurring on a background of systemic and genetic influences. The effect of these influences can be additive, synergistic, or to oppose those of functional adaptation. The study that we report here describes alterations in bone mass and architecture resulting from aging,

and left tibiae from ovariectomized mice, *n* = 6. Dots/bars above the graphs indicate sites of significant difference, *p <* 0.05 following Bonferroni correction. Interactions are the site by ovariectomy interactions calculated by mixed models; ns not significant, \*\*\**p <* 0.001.

changing the mechanical environment, or ovariectomy. It shows that analysis of cortical bone at multiple sites provides more meaningful information than that obtained from analysis of candidate sites because several parameters are altered in a sitespecific manner and so conclusions made about changes in one site may not be meaningfully extrapolated to the remainder of the bone functional unit. This is particularly true in relation to the moments of inertia and curvature, which can only be appreciated when the entire bone is analyzed. These analyses have shown that site specificity analysis provides additional information over conventional approaches and have demonstrated that several physiological contexts, not only mechanical loading, can induce significant site-specific changes in cortical bone mass and architecture. Site specificity analysis quantifies standard cortical parameters precisely and reproducibly. It provides measures, which are directly comparable to those obtained using conventional µCT analysis. Furthermore, using mixed model statistical approaches, it allows comparisons at numerous sites along the bone's length with adequate statistical power.

Most published studies investigating the effects of altered mechanical loading, systemic or genetic interventions on bone only investigate differences at pre-determined sites along the bone's length. This may begin to explain discrepancies in reported responses. For example, the present study using a relatively large sample of 15 mice suggests that axial loading of the mouse tibia tends to reduce marrow area if the 50% site is analyzed, but does not significantly alter marrow area at the 37% site in the same mice. This is consistent with previous studies analyzing the 50% site concluding that axial loading alters the endosteal surface (26–28), whereas studies focused on the 37% site observed no such change (29, 30, 32). It is now widely accepted that changes in local mechanical loading have site-specific effects and it is common practice to compare responsive and less responsive sites along the bone's length (6, 20, 25, 33). The novel methodology described here eliminates the intrinsic bias in selecting candidate sites irrespective of potential inter-group differences in properties such as curvature, and facilitates comparisons by allowing automated determination of cross-sectional parameters at each 1% of the bone's length. Curvature and eccentricity are likely to be influenced by poorly understood interactions between the mechanical, genetic, and hormonal contexts, which ultimately determine bone architecture. For example, in humans, the buckling ratio of the femoral neck is a macro-architectural feature, which deteriorates with age (34, 35) and is influenced by genetic polymorphisms including variants in the paracrine signaling molecule Wnt16 (36–38). Understanding how the magnitude and direction of mechanical cues informs overall bone structure is likely to require integration of biomechanical and mechanistic studies. However, characterization of site-specific and macro-architectural bone phenotypes in transgenic mouse models may begin to elucidate the genetic determinants of bone architecture beyond simply quantifying bone mass at individual sites.

To this end, a few groups have developed systems to map differences in bone mass between wild-type and transgenic mice to specific locations (33, 39, 40). In addition, commercially available software is also able to analyze serial 2D µCT images, as we have previously reported (21). A key advantage of site specificity over currently available software is the ability to automatically analyze images based on percentages of bone length (represented by the absolute number of cross-sectional images for that bone)

rather than analyzing specific lengths. The percentage-based approach accounts for differences in bone length and thereby facilitates alignment between different mice. Another advantage of site specificity is the free availability of the program code provided with this publication, such that it can be adapted to novel approaches. In addition, automated site specificity analysis is very rapid: in our hands, an entire bone can be analyzed on a laptop computer in 4 min (BoneJ analyses being additional) as compared with analysis of a single 100-line section using commercially available software taking 5–6 min. Irrespective of the software used, the general method of analyzing µCT data at multiple sites may begin to inform the mechanisms by which bones achieve and maintain their geometry.

Biological regulation of macro-architectural features is suggested by consistency in properties such as the position of the tibia–fibula junction between different mice. The genetically determined baseline toward which the diaphyses of long bones are programed to develop in the absence of mechanical stimulation is effectively cylindrical (41, 42). Although cylindrical bones would be ideally suited to withstand loading applied in a predictably axial direction, this design is inefficient when the direction of loading cannot be reliably predicted as in the case of footfall. Consequently, bone (re)models toward a compromise between energetically unfavorable increases in mass and functionally integral resistance to fracture. This compromise has been suggested to explain the tendency of bones to be curved; although curvature is inefficient in that it increases strains experienced, it may be beneficial by providing predictability of strain direction (15, 43).

The biomechanical significance of changes in curvature observed following disuse, ovariectomy, and aging, but not following loading in the accustomed axial direction, is poorly understood. As curvature of the mouse tibiae provides a damping effect (15), increased curvature may be an adaptive process in situations of sub-optimal bone mass. This hypothesis is consistent with the observation that the proximal region of the mouse tibia, which has the greatest deviation from the central axis, also has greater mass and is less cylindrical than the distal tibia. The greater polar moment of inertia in the proximal than distal tibia suggests greater resistance to torsional forces (31). Similar regional differences in polar moment of inertia have been reported in the human tibia (44). Finite element modeling may be able to further clarify whether proximal tibial architecture increases distal strain predictability and thereby favors a more cylindrical shape. Given simulated periosteal or endocortical expansion has a greater absolute effect on the calculated polar moment of inertia in the proximal than distal tibia, it is possible that the apparent "unresponsiveness" of polar moments of inertia in the distal tibia following any of the interventions tested may be due to these changes being too small to detect as statistically significant. When effect sizes are small, performing multiple comparisons followed by a Bonferroni *post hoc* test results in Type II errors (at a rate of 1.2% in our simulated 10% change in Ct.Ar with *n* = 6). Future studies aimed at detecting smaller differences may require an increase in sample size to maintain their statistical power. Alternatively less stringent *post hoc* tests (e.g., Tukey's HSD) could be used while accepting a greater proportion of false positive results. Another approach would be to analyze a smaller subset of the data; for example, by only statistically comparing each 5% site to reduce multiple comparisons. Whatever approach is used, graphical representation of site specificity outputs provides a visual representation of the magnitude of differences in different sites, permitting more targeted investigation of sites of potential interest; e.g., by histomorphometry.

Notwithstanding, the distal tibia below the tibia–fibula junction is not "architecturally inert." In a recent study, electrical stimulation of the rat common peroneal nerve resulting in muscle contraction was found to cause bone formation in this distal region of the tibia, correlating with the highest localized strain levels predicted by finite element modeling (45). In the present study, significant differences were observed in distal tibial eccentricity following ovariectomy and aging. Although 2 weeks of disuse has previously been reported to decrease bone eccentricity in growing mice (4), 3 weeks of disuse did not result in significant changes in eccentricity in skeletally mature mice in the present study. Also in growing mice, 10 weeks of ovariectomy significantly increased eccentricity in the distal tibia, potentially reflecting altered loading due to impoverished architecture proximal to the tibia–fibula junction. Aging had the opposite effect, decreasing eccentricity as the distal tibial cortex thinned into a more uniformly circular shape with age.

Aging is associated with reduced bone strength, thinning of the load-bearing cortices, and impoverished bone architecture in mice (19, 26) as in humans (12, 13, 34, 46). Why it is that functionally adapted bones undergo seemingly maladaptive age-related structural changes remains incompletely understood. At the cellular level, age-related deficiencies in the responses of osteoblasts to mechanical strain have been identified (19). At the tissue level, we describe in this study the site-specific reduction in cortical area and increase in medullary area as well as uniform reduction in cortical thickness in the mouse tibia, a pattern of change, which is similar to that seen following neurectomy. These similarities between the effects of aging and neurectomy-induced disuse are consistent with our group's hypothesis that aging reflects a failure of bone's adaptation to mechanical loading such that (re)modeling occurs as in a state of disuse (10).

In summary, automated site specificity analysis of bone structure at multiple sites along a bone's length provides the opportunity to study local and systemic influences on site-specific responses, which alter bone architecture. Application of this system of analysis to the mouse tibia demonstrates that increased

### **References**


or decreased loading, ovariectomy, and aging all cause sitespecific changes in bone, such that the structural alterations achieved by these changes, be they strategic or maladaptive, remain preferentially targeted to specific locations. Even when the imposed (re)modeling stimulus is apparently systemic, as in the case of ovariectomy, alteration in bone structural properties including polar moment of inertia remains site specific. Elucidating the interactions between systemic and mechano-adaptive (re)modeling stimuli, which determine bone architecture, and therefore fracture resistance is expected to be facilitated by integrated analysis of structural responses along the length of the bone functional unit. Site specificity analysis may help refine experimental design by selectively targeting optimally responsive or unresponsive cortical regions. It is anticipated that application of site specificity analysis will begin to identify mechanisms by which genetic or systemic cues enhance or suppresses the local, sitespecific mechanisms underlying the mechanostat to determine the overall mechanical suitability of the bone functional unit.

### **Acknowledgments**

LM is a recipient of an Integrated Training Fellowships for Veterinarians from the Wellcome Trust. GG and PD were supported by the Wellcome Trust when this work was undertaken.

### **Supplementary Material**

The Supplementary Material for this article can be found online at http://www.frontiersin.org/article/10.3389/fendo.2015.00052/ abstract


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Galea, Hannuna, Meakin, Delisser, Lanyon and Price. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Early trabecular development in human vertebrae: overproduction, constructive regression, and refinement**

*Frank Acquaah 1,2, Katharine A. Robson Brown <sup>3</sup> \*, Farah Ahmed <sup>4</sup> , Nathan Jeffery <sup>5</sup> and Richard L. Abel 1,4*

*<sup>1</sup> MSk Laboratory, Department of Surgery and Cancer, Faculty of Medicine, Imperial College London, London, UK, <sup>2</sup> School of Medicine, King's College London, London, UK, <sup>3</sup> Department of Archaeology and Anthropology, University of Bristol, Bristol, UK, <sup>4</sup> Department of Mineralogy, The Natural History Museum, London, UK, <sup>5</sup> Department of Musculoskeletal Biology, Institute of Ageing and Chronic Disease, University of Liverpool, Liverpool, UK*

#### *Edited by:*

*Phil Salmon, Bruker-microCT, Belgium*

#### *Reviewed by:*

*Claire Elizabeth Clarkin, University of Southampton, UK Paula H. Stern, Northwestern University Feinberg School of Medicine, USA*

#### *\*Correspondence:*

*Katharine A. Robson Brown, Department of Archaeology and Anthropology, University of Bristol, 43 Woodland Road, Bristol BS8 1UU, UK kate.robson-brown@bristol.ac.uk*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

> *Received: 30 January 2015 Paper pending published: 04 March 2015 Accepted: 14 April 2015 Published: 01 May 2015*

#### *Citation:*

*Acquaah F, Robson Brown KA, Ahmed F, Jeffery N and Abel RL (2015) Early trabecular development in human vertebrae: overproduction, constructive regression, and refinement. Front. Endocrinol. 6:67. doi: 10.3389/fendo.2015.00067* Early bone development may have a significant impact upon bone health in adulthood. Bone mineral density (BMD) and bone mass are important determinants of adult bone strength. However, several studies have shown that BMD and bone mass decrease after birth. If early development is important for strength, why does this reduction occur? To investigate this, more data characterizing gestational, infant, and childhood bone development are needed in order to compare with adults. The aim of this study is to document early vertebral trabecular bone development, a key fragility fracture site, and infer whether this period is important for adult bone mass and structure. A series of 120 vertebrae aged between 6 months gestation and 2.5 years were visualized using microcomputed tomography. Spherical volumes of interest were defined, thresholded, and measured using 3D bone analysis software (BoneJ, Quant3D). The findings showed that gestation was characterized by increasing bone volume fraction whilst infancy was defined by significant bone loss (*≈*2/3rds) and the appearance of a highly anisotropic trabecular structure with a predominantly inferior–superior direction. Childhood development progressed via selective thickening of some trabeculae and the loss of others; maintaining bone volume whilst creating a more anisotropic structure. Overall, the pattern of vertebral development is one of gestational overproduction followed by infant "sculpting" of bone tissue during the first year of life (perhaps in order to regulate mineral homeostasis or to adapt to loading environment) and then subsequent refinement during early childhood. Comparison of early bone developmental data in this study with adult bone volume values taken from the literature shows that the loss in bone mass that occurs during the first year of life is never fully recovered. Early development could therefore be important for developing bone strength, but through structural changes in trabecular microarchitecture rather than bone mass.

#### **Keywords: vertebra, trabecular, bone growth, development, ontogeny, microcomputed tomography**

### **Introduction**

Gestation, infancy, and childhood are periods of rapid bone growth (increase in size) and development (change in shape). It is generally understood that these early periods of ontogeny may have a significant influence on bone strength and fracture risk during old age (1–4). Intuitively this makes sense because in old age, particularly after the menopause or andropause, bone tissue wastes away, but individuals who develop stronger bone before adulthood will withstand greater age-related loss of bone before a fragility fracture occurs. As yet though very little data have been published regarding development of the material and structural properties that contribute to strength. Further, the small amount of data that have been published to date provides conflicting evidence.

Material properties (i.e., bone mineral density, BMD) have received the most attention and several studies have emphasized the importance of early BMD for adult bone strength. For example, body mass at birth (5) and 1 year (6) was correlated with postmenopausal bone mineral content, which suggests that gestational and infant development may contribute to adult bone strength. However, BMD of long bones actually decreases by about 30% during the first months of life (7) and peak bone mass is not acquired until young adulthood (1–3). Hence, development of BMD after childhood, and even after skeletal maturity, might have more of an impact on bone health in old age.

Mass also contributes to strength, but the development of bone mass has received only moderate attention, in part because of the difficulty of obtaining and analyzing 3D models of bone [e.g., from microcomputed tomography (micro-CT) scans]. The most comprehensive studies of prenatal and infant development have examined the limb bones. During the gestational period (4–9 months), the bone volume fraction (BVF) in the humerus and femur remained constant (8) or increased (9, 10). In infancy, the femur (11) and tibia (12) exhibited a decrease in volume fraction between 0 and 2 years. Therefore, in line with the reports on BMD, the adult pattern bone mass probably develops only after infancy.

Strength is also dependent upon the structure of bone. It is not clear how the development of trabecular microarchitecture contributes to changes in bone mass because these studies used a variety of methods to visualize and measure the trabeculae, making the data difficult to compare. Glorieux et al. (9) and Salle et al. (10) both employed conventional 2D histological histomorphometric methods, sectioning femora under a microscope. Both studies reported that the thickness and number of trabeculae increased during gestation. In contrast, Reissis and Abel (8) applied digital histomorphometric techniques using relatively low-resolution (100 μm) micro-CT scans of femora. Following the conventional studies, trabecular thickness was reported to increase, but number decreased. Postnatally, Ryan and Krovitz (11) and Gosman and Ketcham (12) both analyzed micro-CT scans in 3D to analyze the femur and tibia, respectively. Both limb bones exhibited a decrease in trabecular number between 0 and 2 years, but the femur tended toward decreasing thickness whilst the tibia tended toward an increase. The divergent growth patterns of trabecular number in the femur and tibia could have a huge effect on strength because connectivity is an important determinant of bone strength.

Unfortunately, these studies of mass and microarchitecture provide a fragmented insight into trabecular development as changes in morphology have only been quantified for restricted periods of development. However, the published studies all suggest that human trabecular growth and development is both dynamic and sequential (12, 13–15). This appears to be the case for other species (16). Specifically, the studies published to date and described here indicate that parturition in particular is strongly associated with a reverse in growth trajectory from increasing to decreasing bone BMD and mass, and from a more connected to a less connected structure. The downstream effect of this reduction in bone mass in adulthood/old age is not clear because infant and adult trabecular structure has not been compared.

In order to address the question as to why BMD and mass decrease after birth, if early development is so important for bone strength, more data are required on early bone development, and in particular before and after parturition. Analysis of early growth and development in vertebral bone, a key fragility fracture site in diseases such as osteoporosis, may reveal just how much infancy contributes to adult bone mass. If it is the case that this is an important developmental period in regard to future adult bone health then perhaps it might be possible in future to identify children who are at risk of fragility related skeletal disorders in old age. The aim of this study is therefore to document early bone development from gestation and infancy through to childhood (2.5 years), and then consider whether these periods are important for building adult bone strength.

## **Materials and Methods**

The ontogeny of vertebral trabecular architecture, a key fracture site, was analyzed for the developmental period between 6 months prenatal and 2.5 years postnatal. Microarchitecture was visualized using micro-CT and measured using BoneJ (17), a plugin for ImageJ (National Institutes of Health, USA; http://rsb.info.nih. gov/ij/index.html) and Quant 3D (University of Texas).

### **Sample and Provenance**

A total of 120 vertebrae from complete modern human vertebral columns (C1 to L5) were included in the study. These vertebral series were derived from a 19th century collection of juvenile skeletons of documented age at death held at the Royal College of Surgeons, London. The sample was selected to include three developmental periods: (i) gestational, 6 months – term; (ii) infant, term – 1.2 years; (iii) early childhood, 1.2–2.5 years. The validity of growth curves collected from past populations (including the specimens in this study) has been discussed elsewhere (8).

### **Digital Histomorphometry**

Vertebral columns were scanned at 28 μm<sup>3</sup> voxel size using a Nikon X-Tek HMX-ST 225 system (18). Homologous bone regions were extracted from individual vertebrae by aligning the bodies into a plane defined by three type-II anatomical landmarks (19): the left, right, and anterior apices of the superior cortical rim on the vertebral body. Alignment was followed by extraction of the largest possible spherical volumes of interest (VOI) about the centroid of the vertebral bodies. Trabecular tissue was segmented using a binary global threshold (20) based on the frequency distribution of the gray values (CT numbers) in the scan (21– 23). The threshold value was placed at the trough between the characteristic peaks for air and bone.

Architectural measures were collected including trabecular: BVF (BV/TV); thickness (Tb.Th); number (Tb.N); connectivity density (Conn.D); structure model index (SMI), and degree of anisotropy (DA). The methods and algorithms are detailed elsewhere (17, 21, 24) but a general description follows here. BVF is the amount of bone per unit volume and was calculated by dividing the number of threshold white voxels by the total number of voxels in the VOI (17). Trabecular thickness (Tb.Th) was calculated by fitting the largest sphere possible within a trabecula. The diameter of this sphere represents the thickness of the trabecular element (17). Trabecular number (Tb.N) was measured using the mean intercept length (MIL) method and involved calculating the sum of intersections between trabecular bone and a matrix comprised of a grid of lines (11, 24). Similarly, connectivity density measured the number of connections between trabeculae per cubic millimeter (21). Connectivity density was calculated by dividing the Euler characteristic by the volume of the sphere (17). Structure model index (SMI) measured the proportion of "rod" and "plate"-shaped trabeculae. SMI values ranges from 0 for a structure made of "plate"-shaped trabeculae to 3 representing a primarily "rod"-shaped one. Structures with a mixture of plate and rod shapes are assigned a value between 0 and 3. Anisotropy measured the direction and size of the preferred orientation of trabeculae, calculated as a ratio between the maximum and minimum radii of the MIL ellipsoid (17, 25–27). A value of one reflects a completely anisotropic structure whilst zero reflects an isotropic one. Anisotropic trabecular bone is highly orientated along one direction and is typically highly organized. In comparison, isotropic trabecular bone is not highly orientated and may have a disorganized or randomly orientated structure (12, 28).

### **Statistical Analysis**

Structural growth and development of vertebral trabecular bone was analyzed using a one-way ANOVA with Tukey's *post hoc* carried out using PAST version 2.15 (29). Measurements from each vertebral column (24 bones) were combined into one value for each individual.

### **Results**

Variation in cancellous bone architecture with spinal level is presented in **Figures 1**–**4**. Lumbar vertebrae exhibit significantly lower BVF in relation to cervical at 7 months, term, 1.2 and 2.5 years (one-way ANOVA *p* = 0.038; **Figure 1**). Cervical and thoracic vertebrae exhibit significantly thicker trabeculae than lumbar at 7 months and term (one-way ANOVA *p* = 0.050; **Figure 2**). Cervical and thoracic vertebrae exhibit a significantly higher DA than lumbar at 6 months, 7 months and term, and 2.5 years (one-way ANOVA *p* = 0.001; **Figure 3**). Finally, cervical vertebrae exhibit significantly more numerous trabeculae than lumbar at 6 months, 7 months, 1.2 years, and 2.5 years (one-way ANOVA *p* = 0.008; **Figure 4**).

While the maturational schedule of the regions of the spine appears to vary, the cervical, thoracic, and lumbar vertebrae all exhibited the same general pattern of change in trabecular microarchitecture (**Table 1**). Grouping the bones together, BVF increased prenatally (**Figure 5A**). The rise in BV/TV was significant between 6 and 7 months (*p <* 0.01) but not 7 and 9 months (*p >* 0.05). Increased volume was primarily attributed to growing trabecular thickness between 6 and 7 months (**Figure 5B**) followed by rising number between 7 and 9 months (**Figure 5C**). Elements grew thicker between 6 and 7 months (*p <* 0.01) then thinner from 7 to 9 months (*p <* 0.01). Whilst number decreased between 6 and 7 months (*p >* 0.05) then increased from 7 to 9 months (*p <* 0.01). Connectivity density increased between both 6 and 7 months (*p >* 0.05) and 7 and 9 months (*p <* 0.01; **Figure 5D**). In addition, there was a statistically significant shift from a rodlike to a plate-like structure between 6 and 7 months (*p <* 0.01) and between 7 and 9 months (*p <* 0.01; **Figure 5E**). Anisotropy

decreased during gestational development, between both 6 and 7 months (*p <* 0.01) and 7 and 9 months (*p >* 0.05; **Figure 5F**).

In infancy, there was a significant (2/3rds) decrease in BVF (*p <* 0.01, **Figure 5A**) via decreased trabecular thickness (*p <* 0.01, **Figure 5B**) and number (*p <* 0.01, **Figure 5C**). Connectivity density decreased (*p <* 0.01, **Figure 5D**) and there was a transition from plate-like to rod-like trabeculae (*p <* 0.01, **Figure 5E**), as well as a more highly orientated structure (*p <* 0.01, **Figure 5F**).

During early childhood, BVF remained constant (**Figure 5A**). Trabecular thickness increased (*p <* 0.01, **Figure 5B**) but number decreased (*p >* 0.05, **Figure 5C**). Connectivity density decreased (*p <* 0.01, **Figure 5D**). Trabeculae became more plate-like (*p <* 0.01, **Figure 5E**) and continued to become more highly orientated (*p >* 0.05, **Figure 5F**).

### **Discussion**

The ontogeny of early vertebral trabecular architecture during gestation, infancy, and early childhood was analyzed by spinal level, and then using a cross-sectional approach. The maturational schedule of the regions of the spine vary, with lumbar lagging slightly behind the thoracic region, but a similar pattern

6 months, 7 months, 1.2 years, and 2.5 years (one-way ANOVA *p* = 0.008). Values are from a single sample.


of microarchitectural change is common to all regions. Each cross-sectional developmental stage exhibited distinct patterns of microarchitectural growth and development. Before birth trabecular, BV/TV increased but became relatively isotropic. This pattern was reversed during infancy when BV/TV decreased but the structure became more anisotropic. In early childhood, a new pattern of development appeared, BV/TV remained constant but the DA continued to increase. Both the birth and the end of infancy may be associated with a marked change in growth trajectory of BV/TV.

### **Patterns of Bone Modeling During Early Development**

Variation in the maturational schedule of vertebrae with spinal level was observed; this provides some support for the suggestion that ossification centers first appear in the lower thoracic and upper lumbar regions, and then spread cranio-caudally (30). A similar pattern of microarchitectural change, however, is common to all regions (**Figure 6**). During gestation (6–9 months), there was a steady increase in BV/TV, peaking at term (**Figure 5A**). Term vertebrae actually possessed a higher BV/TV than older infant and childhood specimens. In contrast, during infancy (0–1 years), vertebrae exhibited a significant decrease in BV/TV, which dropped by over 2/3rds from 0.310 at term to 0.079 at 1.2 years (**Figure 5A**). Throughout early childhood (1.2–2.5 years), BV/TV remained constant (**Figure 5A**).

The observed developmental changes in BV/TV are due to variation in the size and number of trabecular elements. Increasing BV/TV during gestation was solely attributable to increased Tb.N (**Figure 5C**) because Tb.Th decreased only slightly (**Figure 5B**). The loss of BV/TV during infancy was due to a reduction in both Tb.Th (**Figure 5A**) and Tb.N (**Figure 5B**). During childhood, BV/TV remained constant because Tb.Th increased (**Figure 5B**) but Tb.N decreased (**Figure 5C**). Hence, gestation is characterized by increasing Tb.N, whereas infancy and childhood is characterized by a decrease in Tb.N. Further, gestation and infancy are characterized by decreasing Tb.Th whilst childhood is characterized by an increase in Tb.Th.

Developmental changes in trabecular microarchitecture suggest that different bone modeling processes occur during gestation, infancy, and childhood. During gestation, the number of nodes increases (Conn.D **Figure 5D**), as might be expected, given the increase in Tb.N. The trabeculae became more plate-like (SMI; **Figure 5E**) suggesting the pores between trabeculae were filled in. Hence, the gestational period seems to be characterized by deposition of bone tissue, the formation of new trabeculae and

**FIGURE 5 | Age-related variation in vertebral trabecular architecture**. Asterix denotes significant differences vs. previous age group (\**p <* 0.05; \*\**p <* 0.01). During gestation (6 months-term) there was a steady increase in BV/TV, peaking at term **(A)**. Increasing gestational BV/TV was attributable to increased Tb.N **(C)** which is also reflected by the increased Conn.D during this time period **(D)**. During infancy (term-1.2 years) vertebrae exhibited a significant decrease in BV/TV **(A)** due to a reduction in Tb.Th **(B)** and Tb.N **(C)**, causing trabecula to become more rod like [SMI **(E)**] and more anisotropic with a dominant inferior-superior axis [DA **(F)**]. During childhood (1.2–2.5 years) there is an increase in BV/TV **(A)** and trabecula return to a more, plate like shape [SMI **(E)**] while structure continues to become more anisotropic, organized along the inferior-superior axis [DA **(F)**].

**FIGURE 6 | Micro-CT reconstructions showing the development of vertebral trabecular microarchitecture of the first lumbar bone during the gestational, infant, and childhood periods**. The size of the bones has been scaled to demonstrate the variation in bone volume fraction.

nodes, as well as the conversion of existing rod-like trabeculae to plates. In contrast during infancy, the number of trabecular nodes decrease (Conn.D; **Figure 5D**) as a result of decreasing Tb.N. The trabeculae became more rod-like indicating the removal of tissue between trabeculae (SMI; **Figure 5E**). Hence, trabecular development in the first year is probably defined by resorption of bone tissue, including entire trabeculae and the perforation of plates. Childhood was the only period that was not characterized by a change in trabecular mass because some trabeculae became thicker (**Figure 5B**) whilst others were removed entirely (**Figure 5C**). Loss of trabeculae probably caused the apparent reduction in Conn.D (**Figure 5D**) whilst deposition of tissue may have affected the increase in SMI (**Figure 5E**) by filling the pores between trabeculae. Thus, childhood development appears to be characterized by selective deposition and resorption of tissue at different sites. Overall then gestational development is characterized by bone deposition, which reverses to resorption during infancy and then gives way to a more balanced pattern of selective deposition and resorption in childhood.

### **Mechanisms Controlling Early Bone Development**

The reversal in growth trajectory (switch from increasing to decreasing BV/TV) after birth could be attributed to a change in the epigenetic mechanisms that control bone growth and development, as well as the loading environment. Uterine bone growth is generally thought to be developmentally preprogramed (31, 32) with some adaptive response to loading (33). However, it is generally assumed that after birth, bone shape is largely influenced by loading (34, 35), probably because parturition marks the beginning of an important functional transition from uterine punching and kicking in the fetal position to habitual weight bearing via sitting, crawling, standing, and walking.

Gestational movements start at about 11 weeks, beginning with sporadic and uncoordinated muscle contractions such as punching and kicking (36–38). More complex movements such as flexion of the limb joints or putting the thumb in the mouth begin at about 5 months and a regular schedule of movement by 6 months of age (39). Yet, vertebral trabeculae do not become more organized (**Figure 5F**) with strongly orientated trabeculae (**Figure 6**), as bone often does when adapting to loads [Sensu (35, 40)]. Trabecular tissue became more isotropic during fetal development (**Figure 5F**). Further, toward the end of gestation, the frequency of movements decreases due to lack of space (41), yet the BV/TV continues to increase toward term. Hence, it seems reasonable to infer that gestational development is programed rather than a response to loading. It has been argued that kicking and punching against the uterine wall (or collisions with other limbs) may represent an intrauterine form of resistance training (7, 42). After birth, the frequency of limb movements increases (43). Yet, the infant's movements typically occur without much resistance, perhaps putting much smaller loads on the skeleton (44). Thus, the loss of bone observed after birth could be argued to occur, at least in part, because bone becomes more responsive to what are actually reduced mechanical loads. However, this is probably not the case because infant bone resorption resulted in a more anisotropic structure (**Figure 5F**) with a dominant inferior–superior orientation (axis from bottom to top in the vertical plane; **Figure 6**). Strongly orientated trabeculae are often referred to as arcades (45) or bundles (46–48), which are typically orientated along the principle strain axis. Ryan and Krovitz (11) reported that the primary arcades in the human proximal femur develop during infancy, with an adult-like pattern appearing at about 2 years. Abel and Macho (46) found that the trabecular bundles in the ilium developed during infancy (*<*6 months) with the adult-like pattern appearing after 12 years. The authors suggested that the appearance of the main trabecular arcades was an adaptation to developing loads as a result of the ongoing emergence and maturation of gait. The loss of tissue during infancy may be essential for developing a highly orientated structure that can resist loads efficiently with minimal bone mass. This would make sense, as it is easier to take away surplus material than add new tissue, which also provides greater phenotypic plasticity.

The continued development of a highly orientated trabeculae structure in the vertebrae during childhood may be a response to developing postural and locomotor loads (49, 50). The dominant loading direction exerted on the vertebrae by sitting, crawling, bipedal standing, and walking were probably along the inferior–superior axis. This is in the form of muscle contraction (49), weight bearing, and ground reaction forces which travel up and down the vertebral column. The structure adapts (34) in response to an altered loading regime; allowing it to become more efficient at resisting compressive loads along an inferior–superior axis. Additionally, the apparent loss of bone that occurred during the first year of life may have served to supply the growing infant with calcium. During pregnancy, maternal calcium levels are depleted leading to reduced concentrations in breast milk that cannot sustain infant levels (51). The observed loss of vertebral bone during infancy coincides with elevated parathyroid hormone levels in comparison to childhood and adolescence (44, 52, 53). Parathyroid hormone increases bone resorption and release of calcium for mineral homeostasis. Therefore, a developmentally programed increase in BV/TV before birth might serve to create a calcium reservoir for the growing infant. Whilst the change from a rod-like to plate-like structure (**Figure 5E**) that occurred during gestation would have created a larger surface area more suited to osteoclastic resorption for mineral release.

### **Term Bone Mass is Never Recovered**

The result of childhood bone remodeling was that trabecular structure continued to become more strongly orientated (**Figure 5F**) along the longitudinal axis (**Figure 6**). During this time, Tb.Th, DA, SMI, and Conn.D metrics all return toward gestational values however the BV/TV and Tb.N never recover. At 2.5 years, both BV/TV and Tb.N metrics are over 2/3rds lower than at term. Comparison with published adult bone volume data indicate that the bone lost during infancy is never completely replaced (**Table 2**). At term, the BV/TV is about 0.310 but



Acquaah et al. Early trabecular development in human vertebrae

Hildebrand et al. (54) reported that in a modern population of adults (22–94 years), the BV/TV was 0.040–0.226. Several other studies investigating contemporary older adults (50 + years), all reported that the highest BV/TV value was between 0.130 and 0.180 (**Table 2**). At term, BV/TV is between 27 and 87% higher than in young or old adults, respectively.

Generally, it is thought that early development is important for bone health, but our data show that the infant period is most important, more so in terms of structure than mass. During this time, bone mass decreases as microarchitecture changes, allowing the bone to adapt to the increased demands of loading during adulthood efficiently. Infant and childhood development could be very important for developing adult-type pattern of trabecular morphology. More developmental data from 2 years to adulthood is required to determine whether this is the case.

### **Conclusion**

Overall the sequential pattern of trabecular growth development appears to be one of gestational overproduction followed by infant "sculpting" of bone tissue during the first year of life, and then subsequent refinement during early childhood. Similar patterns of development have been described in other parts of the body. During embryological hand development, the finger digits are sculpted from the hand paddle by programed cell death occurring in the cells that are no longer needed (62). During brain development, there is an initial overproduction of

### **References**


tissue, which is then refined on the basis of functional activity (63). The data presented in this study appear to show that early trabecular development follows the same basic pattern of overproduction, constructive regression, and refinement. This process probably serves two functional purposes; gestational overproduction creating a calcium reservoir for mineral homeostasis, whilst sculpting of bone during infancy and childhood releases the calcium and allows greater phenotypic plasticity as bone can adapt to the prevailing loading environment. Hence, in line with current thinking (1, 2), early development could be very important for developing bone strength, but through structural changes in trabecular microarchitecture as opposed to mass.

### **Acknowledgments**

Martyn Cooke, Royal College of Surgeons curator for loan of materials. Michael Doube, Royal Veterinary College for assistance with BoneJ. Richard Ketcham, University of Texas for assistance with Quant 3D. Beatriz Pinilla for research assistance in micro-CT scanning.

## **Supplementary Material**

The Supplementary Material for this article can be found online at http://journal.frontiersin.org/article/10.3389/fendo.2015.00067/ abstract

#### **Video 1 | Movie of Figure 6**.

microarchitectural change. *Am J Phys Anthropol* (2009) **138**:318–32. doi:10. 1002/ajpa.20931


*and Informatics*. (Vol. 2), Los Alamitos, CA: IEEE Computer Society (2008). p. 23–7.


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Acquaah, Robson Brown, Ahmed, Jeffery and Abel. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Hypothesis: coupling between resorption and formation in cancellous bone remodeling is a mechanically controlled event**

### *Reinhold G. Erben\**

*Institute of Physiology, Pathophysiology and Biophysics, Department of Biomedical Sciences, University of Veterinary Medicine Vienna, Vienna, Austria*

#### *Edited by:*

*Daniel Chappard, Centre Hospitalier Universitaire Angers, France*

#### *Reviewed by:*

*Jennifer Tickner, University of Western Australia, Australia Roger Brooks, University of Cambridge, UK*

#### *\*Correspondence:*

*Reinhold G. Erben, Department of Biomedical Sciences, Institute of Physiology, Pathophysiology and Biophysics, University of Veterinary Medicine Vienna, Veterinaerplatz 1, Vienna 1210, Austria reinhold.erben@vetmeduni.ac.at*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

> *Received: 27 February 2015 Accepted: 05 May 2015 Published: 20 May 2015*

#### *Citation:*

*Erben RG (2015) Hypothesis: coupling between resorption and formation in cancellous bone remodeling is a mechanically controlled event. Front. Endocrinol. 6:82. doi: 10.3389/fendo.2015.00082* Coupling is the process that links bone resorption to formation in a temporally and spatially coordinated manner within the remodeling cycle. In order to maintain skeletal integrity, it is of crucial importance that the amount of bone resorbed matches the amount of newly formed bone in each remodeling site. Although a number of different explanatory models have been developed, the mechanisms that couple bone resorption and formation in bone remodeling are still a matter of controversy. Here, I propose a model in which coupling is achieved by biomechanical strain sensed by osteocytes within the newly built bone package. In this model, the resorption cavity created by osteoclasts results in mechanical weakening of the structural element, and, thus, in increased strain under constant loading conditions. Subsequent bone formation is initiated by strain-sensitive osteocytes in the underlying bone matrix. After osteoblastic bone formation has started, the newly built osteocyte–osteoblast network detects strain. Once the mechanical strain within the newly built bone structural unit falls below a certain threshold, bone formation stops. In this biomechanical strain-driven model, osteoblasts do not need to "know" how much bone was previously resorbed in a given site. In addition, this model does not require the transfer of any information from bone-resorbing osteoclasts to boneforming osteoblasts, because biomechanical strain "guides" osteoblasts through their job of re-filling the resorption cavity.

#### **Keywords: bone remodeling, coupling, bone formation, bone resorption, disuse**

Bone remodeling is a cyclical renewal process in which, after activation, a quantum of bone is first resorbed by osteoclasts. Thereafter, osteoblasts fill the resorption cavity with new bone in the same place. In contrast, bone formation and bone resorption are not coupled, and occur independently from each other during bone modeling, resulting in resorption or formation drifts which alter bone structure at the microscopic or macroscopic level.

Bone remodeling and bone modeling can both be found in cancellous and cortical bone of higher mammals. In a growing mammalian skeleton, cancellous bone turnover is dominated by modeling, whereas remodeling is the major turnover activity in a mature skeleton (1). In cortical bone, intracortical bone remodeling leaves behind typical microanatomical structures, namely Haversian canals or osteons. In intracortical bone remodeling, osteoclasts and osteoblasts are organized in a complex structure, the so-called basic multicellular unit (BMU). BMUs consist of a cutting cone of osteoclasts, followed by a closing cone lined by osteoblasts, together with connective tissue, blood vessels, and nerves (2). In cancellous bone, it is not entirely clear whether BMUs exist as distinct entities, because the length of the reversal phase appears to be quite variable, at least in postmenopausal osteoporosis (3, 4). The reversal phase is the phase between the end of resorption and the beginning of formation in the remodeling cycle.

It is currently thought that remodeling can be initiated by either stochastic, hormone-driven, or targeted, microdamage-driven, mechanisms. Stochastic remodeling is believed to be under endocrine control, with sex steroids and parathyroid hormone being the main endocrine determinants of bone turnover (5, 6). The purpose of targeted remodeling is to remove microdamage within the bone matrix. However, this distinction between stochastic and targeted remodeling may be arbitrary, because there is currently no proof that both mechanisms operate really independently. In any case, the initial event for initiation of osteoclastic bone resorption in cancellous bone remodeling is likely detachment of bone lining cells from the bone surface, at least in humans (7). Bone lining cells are flat, osteoblast-derived cells covering all quiescent bone surfaces. By detachment of bone lining cells from the bone surface, a canopy is formed under which blood-borne osteoclasts can attach to the bone surface and can start to resorb bone (7). Bone lining cells are able to receive information from osteocytes within the remodeling unit, because they are in contact with underlying osteocytes via gap junctions (8). Osteocytes appear to have a pivotal function not only for detection of microdamage within bone (9) but also for the control of bone turnover via secretion of receptor activator of NFκB ligand (RANKL), an essential cytokine for bone resorption by osteoclasts (10).

The process that links bone resorption to formation in a temporally and spatially coordinated manner within the remodeling cycle is called "coupling." In order to maintain skeletal integrity, it is of crucial importance that the amount of bone resorbed exactly matches the amount of newly formed bone in each remodeling site. A negative bone balance over a longer period of time invariably leads to bone loss and osteoporosis, because a substantial amount of the skeleton is replaced each year in adult humans.

Numerous attempts have been made to explain how the information about the amount of bone resorbed by osteoclasts is transmitted to osteoblasts in the remodeling cycle. It is currently thought that coupling between bone resorption and formation occurs (i) through growth factors stored in the bone matrix, and released during resorption, (ii) through soluble clastokines secreted by osteoclasts, and (iii) through molecules expressed in the cell membrane of osteoclasts [reviewed by Sims and Martin (7)]. Most of our current understanding of the mechanisms involved in coupling comes from experiments in gene-targeted mice. However, mice and rats lack true intracortical, Haversian remodeling (1, 11). Therefore, it is unknown whether there are differences in the coupling mechanisms between intracortical and cancellous bone remodeling.

The current explanatory models of the coupling mechanism are associated with a number of problems. First, none of these models can convincingly explain why the amount of bone formed during the formation phase matches the amount of bone resorbed during the resorption phase. Second, and perhaps more critical is the fact, that in human cancellous bone remodeling, the time span between the end of osteoclastic resorption and the initiation of bone formation is in the range of several weeks (12). Any biochemical signal linking bone resorption to bone formation will have dissipated during this long period of time. Therefore, it is unclear how information is actually transmitted from osteoclasts to osteoblasts. Moreover, a diligently conducted histomorphometric study in human iliac biopsies of patients with postmenopausal osteoporosis revealed a large percentage (~30%) of remodeling cycles that became arrested in the reversal phase (4), suggesting that formation is not always tightly coupled to resorption in cancellous bone remodeling in humans.

These problems led me to hypothesize that coupling in cancellous bone remodeling may simply be a mechanically controlled process within the newly formed bone package. This hypothesis is actually not totally new, because several aspects of it have been described earlier (1, 12–20). However, it is presented here as a synthesis of different elements and in a refined form, taking into account the microanatomy of newly built bone packages. Using finite element models, Huiskes et al. (15) and Smit & Burger (16) provided mathematical descriptions of cancellous bone remodeling and of the potential strain distributions around a resorption cavity in cancellous and cortical bone remodeling, respectively, and suggested that strains sensed by osteocytes within resorption cavities could account for subsequent activation of osteoclasts and osteoblasts. Later on, these theories were further extended by including mathematical models of fluid flow in the osteocyte canalicular system around the resorption tunnel (18), and by simulation models for osteoclast activity (19). These mathematical models may explain why osteoclastic bone resorption proceeds along the loading axis, and why different strain distributions within the resorption cavity may account for spatial differences in the activation of different cell types. However, the latter models did not explicitly address the key feature of the coupling phenomenon, namely that the amount of newly formed bone matches the amount of previously resorbed bone in a given remodeling site. Huiskes et al. (15) proposed that the magnitude of the bone formation stimulus generated by osteocytes located in the underlying bone matrix may determine the number of osteoblasts recruited, and, thus, the amount of bone formed during the formation phase.

The conceptual advance of the current hypothesis is that it provides a plausible and self-regulating mechanism for the control of re-filling of the resorption cavity in cancellous bone remodeling based on the microanatomy of newly formed bone packages. In this model (**Figure 1**), the resorption cavity created by osteoclasts results in mechanical weakening of the structural element, and, thus, in increased strain around the resorption cavity under constant loading conditions. The increased strain is detected by osteocytes in the underlying bone matrix. This part of the hypothesis is supported by finite element models of the strain distribution around resorption cavities (16). When the strain exceeds a certain threshold, the osteocytes initiate subsequent bone formation by secreting osteogenic signals through the canalicular network opened by osteoclastic bone resorption (**Figure 1A**). After osteoblastic bone formation has started, the newly built osteocyte–osteoblast network detects strain, because

the underlying osteocyte canaliculi system is sealed by the cement line (**Figure 1B**). All previous mathematical models have not taken into account that the cement line disrupts osteocyte signaling and also canalicular fluid flow from the underlying bone matrix to the surface. In addition, newly formed bone is less mineralized and has, therefore, different material properties compared with the higher mineralized surrounding old bone. It is likely that the differences in material properties between old and new bone affect strain energy distributions within the newly formed bone package, and, thus, mechanosensing of matrix-embedded osteocytes. Once the mechanical strain within the newly built bone structural unit falls below a certain threshold, bone formation stops. Because wall thickness has to be controlled within a range of a few micrometer to achieve constant trabecular thickness and bone mass, the strain threshold when bone formation stops needs to include the, depending on the species, 3–15 μm wide unmineralized osteoid seam (**Figure 1B**).

This model, which may also be used to further refine the mathematical models of bone remodeling, can explain why bone remodeling restores bone structures more or less in their old shape under unchanged loading conditions, i.e., in a biomechanical steady state. In this biomechanical strain-driven model, osteoblasts don't need to "know" how much bone was previously resorbed in a given site. There is no necessity to transfer any information from bone-resorbing osteoclasts to bone-forming osteoblasts, because biomechanical strain within the newly formed bone package "guides" osteoblasts through their job of re-filling the resorption cavity. Further, this model explains why arrest lines can occur in

bone remodeling units. Remodeling unit-associated arrest lines are frequently found in human (21) and also rat (1) bone sections. Arrest lines are generated when osteoblasts temporarily stop and subsequently resume their bone-forming activity. Based on the proposed model, arrest lines indicate a change in mechanical loading during the formation phase of the remodeling cycle, so that strain transiently falls below the threshold to maintain bone formation.

The proposed model makes a number of predictions which could be used to verify or falsify the model. In line with earlier mathematical models (15, 16, 19), the current model predicts that unloading will result in aborted remodeling and accumulation of resorption cavities. In addition, reduced biomechanical loading would lead to under-filling of resorption cavities. Conversely, increased loading would result in over-filling of resorption cavities. Moreover, shallow resorption cavities may not be filled with new bone by osteoblasts, because the increase in biomechanical strain of the structural element caused by a shallow resorption cavity may not be sufficient to elicit an osteogenic signal by osteocytes. In addition, mechanical disconnection of a structural element by excessive resorption and subsequent complete perforation will cause aborted remodeling, and changes in the material properties (e.g., hypo- or hyper-mineralization) of the newly formed bone will affect wall thickness.

Another interesting aspect of this model is that it could be regarded as a "unifying hypothesis of cancellous bone turnover." It has long been an enigma why cancellous bone modeling and remodeling activities can coexist in a cancellous bone network, and how bone cells differentiate between these two different activities. In agreement with mathematical models reported previously (15, 19), the proposed model suggests that both processes follow the same rules and are just different aspects of the same underlying mechanism. For example, when strain falls below a certain threshold in a given structural element, parts of this element will be removed without subsequent induction of bone formation, resulting in a modeling resorption drift. Re-loading of the same element will induce formation on top of resorption which would be interpreted as remodeling in a histological section. The idea of a strain threshold for initiation of bone formation may also explain the observed lag time between the end of osteoclastic resorption and the initiation of bone formation (12). Especially in individuals with low levels of physical activity, it may take time to maintain strains above the threshold over a certain period of time in a specific remodeling site. The proposed model would predict that the lag time between resorption and formation depends on the biomechanical strain within a given structural element, and should be less in a high strain environment. In addition, it is possible that the thresholds may be modulated by endocrine signals (*vide infra*).

Is there any evidence for the validity of this model? In fact, there is. In a scanning electron microscopic study in lumbar vertebrae of normal subjects of different ages, Mosekilde (13) observed that resorption cavities were not filled with new bone on trabecule which lost 3D connection, i.e., unloaded trabecule. Moreover, partial unloading decreased mineral apposition rate, bone formation rate, and wall thickness in the presence of unchanged osteoclast numbers in a rat hindlimb immobilization model (**Figure 2**).

All these experimental findings are predicted by the proposed model.

The current hypothesis is dealing with the larger perspective of the physiological regulatory mechanisms underlying the coupling mechanism in cancellous bone remodeling, not with the actual biochemical signals involved in cell–cell communication. The nature of these biochemical signals is still controversial. However, it is highly likely that the Wnt pathway plays an important role for the osteocyte- or lining cell-driven initiation of bone formation in the proposed model. For example, treatment of rats and monkeys with anti-sclerostin antibodies induces modeling drifts on quiescent bone surfaces together with overfilling of remodeling units, thus, mimicking the effects of biomechanical over-loading (22). Therefore, the strain-induced downregulation of the Wnt inhibitors sclerostin or dickkopf-1 may be an important part of the physiological signaling mechanism in coupling. Thus, the nature of the "osteogenic" signal could in fact be a downregulation of inhibitory signals. Interestingly, mathematical models assuming that mechanical loading inhibits osteocytes from inhibiting bone formation via sclerostin secretion have been shown to be able to produce a load-aligned trabecular structure (23). It is also conceivable that the osteocyte-derived coupling signals are subject to hormonal modulation. For example, tight coupling was observed in patients with primary hyperparathyroidism (3). In this context, it is known that intermittent parathyroid hormone treatment suppresses sclerostin expression in bone (24). Therefore, it is possible that endocrine signals are able to modulate the osteocytic secretion of coupling factors.

The provocative aspect of the proposed model is that it challenges the notion of hemiosteonal remodeling with predetermined BMUs in cancellous bone remodeling. Similar to strain-driven mathematical models of bone remodeling reported earlier (15, 16, 19), osteoclastic bone resorption and osteoblastic bone formation are not directly associated, but rather indirectly linked through mechanosensing osteocytes in the proposed model.

Whether loading might also influence the re-filling of resorption tunnels in Haversian remodeling in cortical bone is currently unknown. The same principles may apply, and mathematical models have been able to modulate both cancellous and intracortical bone remodeling, using similar parameters (19). However, an important difference between cancellous and cortical bone remodeling is that re-filling of the resorption tunnel in cortical bone remodeling is more or less complete, just leaving the Haversian canal in the center of the osteon. Unlike cancellous bone remodeling, the concentric contraction of the bone-forming surface during the formation phase in Haversian remodeling may simply result in mutual steric inhibition of bone-forming osteoblasts, so that bone formation automatically comes to a halt. Because there is no evidence that osteoblasts can move, most of the osteoblasts at a given position of the closing cone must undergo apoptosis. Therefore, theoretically a self-regulating process controlling the amount of bone formed by osteoblasts in the closing cone would not be necessary in Haversian remodeling. It might be enough just to initiate the process by the release of straindependent osteogenic signals after osteoclastic bone resorption. On the other hand, based on the concepts presented here for cancellous bone, strains sensed within the newly built bone matrix might also be necessary for maintenance of bone formation and complete re-filling of resorption tunnels in Harversian remodeling. This notion is supported by an earlier study in non-human primates, showing that transient unloading inhibited intracortical bone formation (25). Similar to cancellous bone remodeling, the lower degree of mineralization of the newly built bone matrix compared with the surrounding cortical bone might influence mechanosensing of osteocytes within the closing cone. Clearly, stringent experiments need to be designed to prove or disprove the proposed model of coupling in cancellous and possibly cortical bone remodeling.

### **Acknowledgments**

The author thanks all his team members for their contributions to the concepts described in this paper which developed over many years. The work of the author is currently supported by grants from the Austrian Science Fund (FWF, P21904-B11, P24186-B21, I764-B13, I734-B13).

### **References**


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Erben. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# **Using micro-CT derived bone microarchitecture to analyze bone stiffness – a case study on osteoporosis rat bone**

*Yuchin Wu<sup>1</sup> , Samer Adeeb1,2 and Michael R. Doschak 1,3 \**

*<sup>1</sup> Department of Biomedical Engineering, University of Alberta, Edmonton, AB, Canada, <sup>2</sup> Department of Civil and Environmental Engineering, University of Alberta, Edmonton, AB, Canada, <sup>3</sup> Faculty of Pharmacy and Pharmaceutical Sciences, University of Alberta, Edmonton, AB, Canada*

#### *Edited by:*

*Phil Salmon, Bruker-microCT, Belgium*

#### *Reviewed by:*

*Phillip Trevor Newton, Karolinska Institutet, Sweden Simone Tassani, Universitat Pompeu Fabra, Spain*

#### *\*Correspondence:*

*Michael R. Doschak, Faculty of Pharmacy and Pharmaceutical Science, University of Alberta, 2-020J, Katz Group Centre for Pharmacy and Health Research, 11361-87 Avenue Edmonton, AB T6E 2E1, Canada mdoschak@ualberta.ca*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

> *Received: 03 February 2015 Paper pending published: 04 March 2015 Accepted: 03 May 2015 Published: 20 May 2015*

#### *Citation:*

*Wu Y, Adeeb S and Doschak MR (2015) Using micro-CT derived bone microarchitecture to analyze bone stiffness – a case study on osteoporosis rat bone. Front. Endocrinol. 6:80. doi: 10.3389/fendo.2015.00080* Micro-computed tomography (Micro-CT) images can be used to quantitatively represent bone geometry through a range of computed attenuation-based parameters. Nonetheless, those parameters remain indirect indices of bone microarchitectural strength and require further computational tools to interpret bone structural stiffness and potential for mechanical failure. Finite element analysis (FEA) can be applied to measure trabecular bone stiffness and potentially predict the location of structural failure in preclinical animal models of osteoporosis, although that procedure from image segmentation of Micro-CT derived bone geometry to FEA is often challenging and computationally expensive, resulting in failure of the model to build. Notably, the selection of resolution and threshold for bone segmentation are key steps that greatly affect computational complexity and validity. In the following study, we evaluated an approach whereby Micro-CT derived grayscale attenuation and segmentation data guided the selection of trabecular bone for analysis by FEA. We further correlated those FEA results to both two- and three-dimensional bone microarchitecture from sham and ovariectomized (OVX) rats (*n* = 10/group). A virtual cylinder of vertebral trabecular bone 40% in length from the caudal side was selected for FEA, because Micro-CT based image analysis indicated the largest differences in microarchitecture between the two groups resided there. Bone stiffness was calculated using FEA and statistically correlated with the three-dimensional values of bone volume/tissue volume, bone mineral density, fractal dimension, trabecular separation, and trabecular bone pattern factor. Our method simplified the process for the assessment of trabecular bone stiffness by FEA from Micro-CT images and highlighted the importance of bone microarchitecture in conferring significantly increased bone quality capable of resisting failure due to increased mechanical loading.

**Keywords: micro-computed tomography, microstructural parameters, finite element analysis, bone stiffness**

## **Introduction**

Micro-computed tomography (Micro-CT) is a common tool in most bone biology laboratories and acknowledged as a valuable technique for investigating the microarchitecture of bone, for both *in vivo* and *ex vivo* applications (1–4). The Micro-CT derived images from reconstructed data accurately represent bone microstructural parameters for quantitative assessments (5). In addition, as the X-ray passes through the sample, the relative linear attenuation of the sample is represented on the image projections as values of gray scale (6). Corresponding X-ray attenuation on reconstructed images can thus be converted to bone mineral density (BMD) by calibrating the values with a tissue "phantom" of known mineral density. For trabecular bone, the structural parameters derived from Micro-CT data are based on traditional static bone histomorphometry, which evaluated the thickness, connectivity, distribution, and spacing of the trabecule. It is important to acknowledge that image segmentation methods can also influence the correlation of structural parameters to the stiffness of the trabecular bone structure being modeled, as described in several key finite element analysis (FEA) publications (7–9). Our key message remains in the added utility of initial selection of an appropriate region of bone for analysis (informed by grayscale attenuation) for building that FEA model. Measured parameters by Micro-CT including bone volume/tissue volume (BV/TV; %), trabecular number (Tb.N; mm-1), trabecular thickness (Tb.Th; mm), and trabecular separation (Tb.Sp; mm) have all been successfully utilized to assess the microarchitecture of rodent trabecular bones (5). In addition, there are other parameters that can also be used to represent bone microarchitectural complexity, such as fractal dimension (FD), trabecular bone pattern factor (Tb.Pf; 1/mm), degree of anisotropy (DA), and connectivity density (Conn.Dn; 1/mm<sup>3</sup> ) (10). Although those parameters are well established with their own definitions, they remain indirect indices and require other computational tools to interpret bone stiffness and potential for mechanical failure (i.e., fracture) (11, 12). Resultantly, in bone research, the numerical method, finite element, was applied to assist in analyzing bone stiffness (13, 14).

Finite element method is founded on numerical methods for calculating the mechanics in various fields in industries. For bone analysis, the procedure starts by acquiring computed tomography images to build the geometry of the bone, which may include both trabecular and cortical components. Each voxel in the image stacks can be converted to an element in the finite element mesh. The material properties required for the analysis can also be inferred from the gray scale values on the images to enable a more accurate representation of the local stiffness of different bone structures (15–17). This process has been successfully implemented in many animal models (18), including modeling the entire rat vertebra, the vertebral body (without the top and bottom growth plates), and the trabecular bone compartment alone (1, 19). Nevertheless, the procedure from image segmentation to finite element model can prove challenging to the non-engineer, and inappropriate bone modeling techniques may result in redundant calculation in FEA. Thus, our aim was to simplify the selection of Micro-CT derived bone geometry for application to FEA and subsequently correlate the FEA results to bone microstructural and density based material parameters as an initial rudimentary assessment of bone stiffness (20). Although previous investigators have plotted the 2D curves of bone structural parameters (21), we show in this manuscript the additional utility of using the curve of 2D structural parameters of bone from Micro-CT image data to inform the selection of bone for subsequent FEA. In addition, using this approach, regional trabecular structural parameters become directly correlated to regional trabecular stiffness from FEA, further informing the impact of the given experimental intervention. In addition, using this approach, regional trabecular structural parameters become directly correlated to regional trabecular stiffness from FEA, further informing the impact of the given experimental intervention.

Our study hypothesis was that structural and density based parameters of trabecular bone microarchitecture from Micro-CT can be correlated to bone stiffness.

## **Methods**

### **Animals**

About 20, 3-month-old female, Sprague-Dawley rats were obtained from Charles River. The protocol pertaining to all procedures and aspects of the study was approved by the University of Alberta animal care and ethics committee. Before arrival, to induce an osteoporosis (OP)-like increase in bone resorption, ten rats were ovariectomized (OP-OVX); the remaining ten rats were sham-operated as normal controls (OP-Sham). All rats were euthanized 12 weeks after OVX and/or sham surgery and lumbar vertebrae 4–6 (L4–L6) were dissected. The muscle and soft tissues connected to vertebrae were removed, and samples were wrapped in paper towels dampened with phosphate buffered saline and stored frozen at -30°C for subsequent Micro-CT scan.

### **Micro-CT Imaging**

The excised vertebrae, L4, L5, and L6, were thawed at room temperature for 2 h before scanning. All vertebrae were imaged using a benchtop Micro-CT imager (SkyScan 1076; Bruker-MicroCT, Kontich, Belgium) at 18 μm voxel image resolution with 70 kV, 100 μA, and a 1.0 mm aluminum filter. Projection images were reconstructed using bundled vendor software (Nrecon 1.6.1.5; Bruker-MicroCT, Kontich, Belgium). Dataviewer 1.4.3 software was employed to orient the cross sectional images parallel to the transaxial plane. A cylindrical region of interest (ROI) was segmented from trabecular bone within the vertebral body. The diameter of each cylinder was constrained by the endocortical bone margin enclosing the vertebra. The vertebral growth plates at each end were used to landmark the top and bottom segmentation boundaries. We utilized CTAn 1.11.6.0 (Bruker-MicroCT, Kontich, Belgium) to morphometrically analyze the trabecular bone cylinder. In order to calibrate BMD with Hounsfield units (HU), two hydroxyapatite [Ca10(PO4)6(OH)2] phantoms (Computerized Imaging Reference Systems Inc., Norfolk, VA, USA) with BMD 0.250 and 0.750 g/cm<sup>3</sup> were included in the calculation based on the assumption that HUair = -1000 and HUwater = 0. The threshold, 70 ~ 255, was chosen to segment out bone on 8-bit (0 ~ 255 gray level) bitmap (BMP) images. The measured structural parameters were BV/TV, BMD, Tb.Th, Tb.Sp, FD, Tb.Pf, DA, and Conn.Dn. Those parameters were calculated three dimensionally (3D) based on the volume of the cylinder. Except for Conn.Dn, all structure parameters were also calculated two dimensionally (2D) from cross sectional images slice by slice, along the longitudinal direction from caudal to cranial (**Figure 1**), and graphed with respect to the normalized height of the cylinder.

**FIGURE 1 | Steps for calculating stiffness of trabecular bone cylinder of rat vertebra**.

The results were compared between groups and different levels of vertebrae in the same group.

### **Finite Element Analysis**

The results of the Micro-CT structural parameters indicated that the largest difference between groups were in the 40% portion of the cylindrical length from the caudal side. Thus, it was selected for FEA (**Figure 1**). The cross sectional images of this region were imported into Mimics 14.11 (Materialise, Leuven, Belgium) to build finite element meshes. Since the voxel of the Micro-CT data was 18 μm (*≈*17.156 μm), that dimension was used to build voxel mesh (hexahedral elements) for FEA to maintain the details of the microarchitecture. To avoid losing or distorting microarchitecture of trabecule, no resizing or smoothing was undertaken prior to meshing. After creating hexahedral elements (hex8), the meshes were exported as an input file (\*.inp) format of ABAQUS (SIMULIA, Providence, RI, USA). The remaining preprocessing to include material properties and boundary conditions were done in ABAQUS/CAE 6.10-1. Since bone is expected to behave elastically and failure expected with small strain (22), the analysis was set as linear and elastic with Young's modulus (E) = 24.5 GPa and Poisson's ratio (ν) = 0.3 (23). All vertebrae were assigned with the same material property, as no significant difference of material properties was measured between normal and OVX bone (23). The bottom of the cylinder was fixed and a displacement (Δd), 5% of the height of the cylinder (strain = 0.05), was applied to the top of the cylinder. The stiffness (k) of the bone was calculated by dividing the summation of the vertical reaction

forces (R) on the fixation to the applied displacement. Since the sizes of individual vertebrae varied, the stiffness (which comprises information of geometry and material property) would also vary with size. Therefore, values of k were further normalized with the theoretical stiffness (k*′* ), which was calculated from the size of each cylinder by assuming that the cylinder was solid (**Figure 2**). The index of bone stiffness was represented as k/k*′* .

### **Statistics**

PASW® Statistics 17.0 was used for statistical evaluations. The twoindependent-samples test, Mann–Whitney *U* test, was used to find differences between individual groups on the same level of vertebra. This method was applied to both the results of Micro-CT structure parameters and the stiffness calculated by the FEA. The difference among three levels (L4, 5, and 6) of the same group was tested by Kruskal–Wallis *H* test. The Mann–Whitney *U* test was done to understand the difference between two individual vertebral levels. Results were shown as mean *±* SD. An asymptotic value (*p*-value) <0.05 was used to evaluate the significance of differences. The correlation between Micro-CT structural parameters and the index of stiffness was described with Pearson correlation coefficients (*r*). The significance level was set at 0.05 (Sig. *<*0.05) with the two tailed test.

### **Results**

### **3D Structural Parameters**

For the same level (i.e., Micro-CT slice) of vertebrae, most computed indices showed significant differences between OP-Sham and OP-OVX, except Tb.Th and DA. No significant difference of Tb.Th was found between OP-Sham and OP-OVX in the vertebrae of this OP-like rat model; however, the formation of new trabecular bone at the vertebral growth plates in these juvenile (3-monthold) rats proceeded as for normal (OP-Sham) rats. Since OP is a systemic metabolic disease, all vertebrae should experience the

#### **TABLE 1 | 3D structure parameters of rat lumbar vertebrae calculated from Micro-CT images**.


*Data expressed as mean and SD.*

*p–value <0.05 means significantly different from other; p–value >0.05 was shaded with gray comparing to the other treatment.*

*<sup>a</sup>Significantly different from L4 under the same treatment.*

*<sup>b</sup>Significantly different from L5 under the same treatment.*

*<sup>c</sup>Significantly different from L6 under the same treatment.*

*BV/TV, bone volume/tissue volume; BMD, bone mineral density; Tb.Th, trabecular thickness; Tb.Sp, trabecular separation; FD, fractal dimension; Tb.Pf, trabecular bone pattern factor; DA, degree of anisotropy; Conn.Dn, connectivity density.*

same changes after OVX. Therefore, in addition to comparing the results of each vertebral body at the same level, the morphometric changes from L4 to L6 in the same group were analyzed to provide further evaluation of OP-like bone change.

The different patterns of structural parameters among the same group and between the two groups (i.e., OVX-Sham vs. OVX-OP) were investigated to evaluate the spatial effect of OP in successive lumbar vertebral body segments (i.e., L4–L6). When comparing values between different levels of vertebrae in OP-Sham, no obvious trends between L4 and L6 were observed. However, in both groups, BMD, Tb.Th, and DA showed significant differences between L4 and L6 (L6 *>* L4), as well as L5 and L6 (L6 *>* L5). FD in OP-OVX showed significant differences among vertebral levels (**Table 1**) with the values increasing with sequential vertebral segment. This phenomenon was not seen in OP-Sham.

### **2D Structural Parameters**

The structural parameters of image slices were calculated and plotted as curves to show the variations along the caudal-cranial direction. As the curves in the same group showed the same shapes, the curves of L4, L5, and L6 were averaged to represent each group. The curves of BV/TV, BMD, and FD were measured to decrease toward the middle; the curves of Tb.Sp and Tb.Pf increased toward the middle (**Figure 3**). In accordance with the observation from 3D models built from Micro-CT images, this reflected the vascular foramina present in the vertebral body. Regarding the subtractions between groups, only Tb.Sp showed the opposite phenomena, since its measurement was magnified by the anatomical structure, namely the vascular foramen (**Figure 4**).

### **Finite Element Analysis**

Based on the reaction force extracted from FEA, the normalized stiffness values (k/k*′* ) were 0.383 *±* 0.092 for OP-Sham and 0.141 *±* 0.053 for OP-OVX. This stiffness index of OP-OVX was significantly less (*p <* 0.05) than that for OP-Sham. The Micro-CT structural parameters of this 40% length of the trabecular cylinder were also calculated in 3D. The results showed that Tb.Th. was not different (*p >* 0.05) between OP-Sham and OP-OVX. The Pearson's correlation coefficients showed that BV/TV, BMD, and FD were positively correlated with the index of stiffness, and the correlation was significant at the 0.05 level. On the contrary, Tb.Sp and Tb.Pf were negatively correlated with the index of stiffness. The values of Conn.Dn were significantly different between groups; yet, they showed no significant correlation to the calculated stiffness (**Table 2**).

### **Discussion**

Our study details an alternative approach for FEA of bone, whereby Micro-CT derived grayscale attenuation and segmentation data serve initially to guide the selection of trabecular bone regions for analysis by FEA. Thus, instead of modeling the entire bone with FEA by default (currently the most common approach employed), our method will determine sample regions with "distinct differences" in the trabecular structure (informed

**FIGURE 3 | 2D structural parameters, BV/TV(%), Tb.Pf(1/mm), FD, Tb.Sp(mm), BMD(g/cm<sup>3</sup> ), of the cylinder trabecular bone along the longitudinal axis from the caudal to cranial direction**. x-axis was normalized by the height of each vertebral body (trabecular cylinder) (x-axis: Caudal to Cranial; 0–1).

by grayscale attenuation and segmentation data) that will be modeled to subsequently identify the relationship between trabecular structure and stiffness. Future implementation of this approach would greatly simplify the assessment of trabecular bone structural parameters and the distribution of stiffness in regional bone locations (e.g., for high-throughput preclinical applications, such as gaging the antiresorptive efficacy of a novel OP drug), without the need for exhaustive computational modeling of the entire bone structure.

This study recapitulates the importance of trabecular bone microarchitecture in conferring structural strength to otherwise similar trabecular bone geometry. In this study, hormonal depletion induced in female rats (secondary to ovariectomy surgery) served to imbalance the bone remodeling process toward systemic bone resorption and decreased bone volume. Despite that systemic bias, no significant difference of Tb.Th between OP-Sham and OP-OVX groups were measured in vertebral trabecular bone using Micro-CT, suggesting that microarchitectural differences were responsible for the significantly reduced structural strength measured in OP-OVX rats.

**TABLE 2 | 3D structural parameters and index of stiffness of the 4/10 length L4**.


*Data expressed as mean and SD.*

*p–value <0.05 means significantly different from the other group; p–value >0.05 was shaded with gray.*

*The significant level of Pearson's correlation coefficient (r) was set at 0.05 and means no significant correlation to k/k′ ; the value (Sig.) >0.05 was shaded with light up diagonal. BV/TV, bone volume/tissue volume; BMD, bone mineral density; Tb.Th, trabecular thickness; Tb.Sp, trabecular separation; FD, fractal dimension; Tb.Pf, trabecular bone pattern factor; DA, degree of anisotropy; Conn.Dn, connectivity density; k/k′ , normalized stiffness value.*

### **Resolution Issues of Images in Micro-CT**

The selection of resolution for Micro-CT imaging is a critical parameter that, if not chosen adequately, would have high impact on the analysis results. The average trabecular thickness of OP-Sham obtained from the analyzed images was 108 *±* 9 μm (**Table 1**). That number was obtained from images scanned with a target resolution of 18 μm. In theory, other resolutions (9 and 35 μm) available through the Micro-CT imager in our laboratory are feasible for scanning rat vertebra. However, the 18 μm resolution meant that each trabecula would contain six voxels along the thickness direction, which is adequate to capture the complex nature of the trabecular structure. In addition, the finite element mesh with six elements per trabecular thickness is considered adequate for the model accuracy. It is worth noting that a possible advantage associated with scanning with a higher resolution is to observe possible microcracks in the trabecular bone. Nevertheless, scanning our samples with a higher resolution (9 μm) was not practical because of the cumbersome file size of the resulting image and the associated doubling of scanning time.

### **Threshold Selection for Segmenting Bone in Image Processing**

Besides the resolution, the threshold for segmenting bone is also essential, since the results of the variation along the trabecular longitudinal axis may change with different thresholds. Thus, consistent threshold should be applied in the same study to exclude differences resulting from variations in thresholds. In general, for calculating the BMD of trabecular bone among different groups, if the HU differences were not significant between groups, the trend of BMD variations should be very similar to those of BV/TV in the 2D curves (**Figure 3**). This was observed from our results and served as evidence of no significant HU variations between normal and OP bone.

### **Bone Selection for FEA and the Advantage of Our Model**

In rodents, histomorphometric measurements of bone remodeling are complicated by concurrent longitudinal bone growth, which expands both the bone diameter as well as trabecular bone at the growth plates. The bone remodeling, which was influenced by ovariectomy, would best be evidenced on trabecular surfaces. Our analysis of segmented trabecular cylinders excluded bone near the growth plates as well as cortical bone. Excluding cortical bone in FEA to evaluate stiffness allowed us to isolate the effect of OP on the stiffness of the trabecular bone. Furthermore, in our model, we analyzed 40% of the trabecular cylinder starting from the caudal side. The values of DA were also reduced by selecting only 40% of the trabecular cylinder. It means the structure

### **References**


changed from being anisotropic toward isotropy (**Tables 1** and **2**) (24, 25). The advantage of using only this selection is in avoiding the foramina effects which may reduce the difference between groups. Furthermore, the model with 40% of the length has fewer elements than the model of entire vertebra or trabecular bone. This led to less computational time for the FEA.

### **Material Properties of Finite Element Model**

A limitation to our study was ignoring the mapping of material properties according to different HU or gray values. The Young's modulus in each trabecula can be correlated to values of HU (26). However, several studies on rat bone of OVX (0.91 *±* 0.13 and 21.01 *±* 2.48 GPa) and sham (0.90 *±* 0.09 and 22.03 *±* 2.44 GPa) groups showed no significant difference in the hardness and elastic modulus (23, 27). Hence, the same material properties were assigned for the two groups in this study. If further mechanical analysis was required, such as permanent deformation, viscosity, or time effect and crack propagation, then material properties would need to be assigned specifically and non-linear FEA could be applied.

### **Conclusion**

In this paper, we presented a simple process by which the stiffness of trabecular bone can be inferred from structural parameters, without the need for computationally exhaustive FEA models. The 2D structural parameters calculated from Micro-CT analysis reflected the difference in the geometry of the bone between the sham and the OP group. The changes in the following Micro-CT parameters (BV/TV, BMD, Tb.Sp, FD, Tb.Pf, and Conn.Dn) between the two groups were examined and compared with the stiffness obtained from FEA. Other than for Conn.Dn, the values of BV/TV, BMD, FD, Tb.Sp, and Tb.Pf were shown to be related to the structural stiffness. This approach can be readily applied to assess trabecular bone structural changes to quickly infer stiffness, particularly during longitudinal "*in vivo*" studies of laboratory rodents, or to assess the effects of drug treatments.

### **Acknowledgments**

This research was funded by the Alberta Osteoarthritis Team Grant from the Alberta Innovates – Health Solutions (AIHS) Interdisciplinary Team Grant program.

with transiliac bone biopsy. *Osteoporos Int* (2010) **21**(2):263–73. doi:10.1007/ s00198-009-0945-7


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Wu, Adeeb and Doschak. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

# Modalities for visualization of cortical bone remodeling: the past, present, and future

*Kimberly D. Harrison and David M. L. Cooper\**

*Department of Anatomy and Cell Biology, University of Saskatchewan, Saskatoon, SK, Canada* 

Bone's ability to respond to load-related phenomena and repair microdamage is achieved through the remodeling process, which renews bone by activating groups of cells known as basic multicellular units (BMUs). The products of BMUs, secondary osteons, have been extensively studied via classic two-dimensional techniques, which have provided a wealth of information on how histomorphology relates to skeletal structure and function. Remodeling is critical in maintaining healthy bone tissue; however, in osteoporotic bone, imbalanced resorption results in increased bone fragility and fracture. With increasing life expectancy, such degenerative bone diseases are a growing concern. The three-dimensional (3D) morphology of BMUs and their correlation to function, however, are not well-characterized and little is known about the specific mechanisms that initiate and regulate their activity within cortical bone. We believe a key limitation has been the lack of 3D information about BMU morphology and activity. Thus, this paper reviews methodologies for 3D investigation of cortical bone remodeling and, specifically, structures associated with BMU activity (resorption spaces) and the structures they create (secondary osteons), spanning from histology to modern *ex vivo* imaging modalities, culminating with the growing potential of *in vivo* imaging. This collection of papers focuses on the theme of "putting the '*why*' back into bone architecture." Remodeling is one of two mechanisms "*how*" bone structure is dynamically modified and thus an improved 3D understanding of this fundamental process is crucial to ultimately understanding the "*why*."

#### Keywords: basic multicellular unit, bone, remodeling, micro-CT, synchrotron

### Introduction

Bone tissue is three-dimensionally (3D) complex in structure and undergoes continual dynamic change. Despite its rigid structure, it is remarkable in its ability to adapt in response to mechanical stimuli associated with loading and to microdamage endured throughout life. Since Clopton Havers' (1) description of "Haversian" canals and iconic works describing microscopic bone structure/function relationships (2, 3), it has been well appreciated that bone renews itself via the turnover of tissue, which we have come to know as "remodeling" (4). Remodeling is critical for maintaining healthy bone tissue; however, it can also lead to age-related bone loss through an imbalance between osteoclastic (bone resorption) and osteoblastic (bone formation) activity. A progressive deficit in bone formation leads to enlarged osteonal canals and thus increased cortical porosity (5–7). Ultimately, this contributes to bone's fragility which is characteristic of osteoporosis (7, 8). Related fractures are significant events in the lives of those afflicted and are frequently associated with serious complications and even

### *Edited by:*

*Phil Salmon, Bruker-microCT, Belgium*

#### *Reviewed by:*

*Paul Zaslansky, Charité Hospital, Germany Kostas Verdelis, University of Pittsburgh, USA*

#### *\*Correspondence:*

 *David M. L. Cooper, Department of Anatomy and Cell Biology, University of Saskatchewan, 107 Wiggins Road, Saskatoon, SK S7N 5E5, Canada dml.cooper@usask.ca*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

> *Received: 03 February 2015 Accepted: 24 July 2015 Published: 11 August 2015*

#### *Citation:*

*Harrison KD and Cooper DML (2015) Modalities for visualization of cortical bone remodeling: the past, present, and future. Front. Endocrinol. 6:122. doi: 10.3389/fendo.2015.00122*

mortality. With increasing life expectancy, osteoporosis and other degenerative diseases of bone are a growing concern for health care systems worldwide (9). As such, study of the spatio-temporal regulation of remodeling is a topic of great significance within bone biology with the potential to impact many lives.

First described by Frost (10), basic multicellular units (BMUs) are the cellular groups responsible for carrying out the remodeling process. In cortical bone, this is achieved through the localized resorption of a cylindrical space (osteoclastic "cutting cone") followed by concentric infilling of new tissue (osteoblastic "closing cone") (**Figure 1**). The resulting structure is referred to as a secondary osteon (synonymous with "*Haversian system*"). As BMUs organizationally lie between the level of the cell and that of the tissue, Frost referred to them as "intermediary" (11). Despite decades of study, our understanding of the intermediary organization of bone remains rudimentary. BMUs are temporary collections of cells brought together to turnover a discrete packet of bone. Their course through bone tissue, their "behavior," is challenging to directly probe, and thus much of our understanding has been inferred from osteon morphology. The orientations of secondary osteons appear to reflect principal stresses (12–16), and thus it has been hypothesized that the progression of BMUs is influenced by mechanical stimuli. This is not surprising as the

FIGURE 1 | Illustration of a BMU showing the classic 'cutting' and 'closing' cone morphology. Osteoclasts located in the BMU's cutting cone are attracted to areas of damaged bone indicated by the microcrack (26) as well as change in the canalicular network (black lacunae indicative of osteocyte apoptosis) caused by mechanical stimuli such as cyclic loading (based on (24)).

two-dimensional (2D) geometry of secondary osteons has been linked to the function of the bones in which they are found (17, 18) and resultant mechanical strains (19). Additional examples include intra-element regional (i.e., anterior, posterior, lateral, medial) variation in osteon morphology (20, 21) and a relation between osteon size and weight observed in humans (22). To explain the link between mechanics and BMU orientation, computational (*in silico*) modeling has looked to stimuli such as localized strain (23) and strain-related fluid flow (24) around cutting cones. Such *in silico* models continue to become increasingly sophisticated, extending into the realm of simulation (25). All models, however, have relied upon highly idealized BMU morphology, and it is unclear how compatible their findings are with the more complex 3D morphologies which have been reported.

Another hypothesis regarding BMU regulation holds that their activities are spatially "targeted" (27) to remove damage manifested as microcracks (16, 28–30). While debate remains over the degree to which remodeling is targeted vs. untargeted (26), if targeting occurs, even for a portion of remodeling events, there must be mechanisms which actively steer BMUs toward damaged areas (**Figure 1**). This has been envisioned as attraction toward an "effective damage removal area," which provides the means by which the osteoclasts of a BMU are drawn toward microdamage (15). Although it has been demonstrated that remodeling-related resorption spaces are associated with microcracks (31), active steering has yet to be empirically demonstrated. It is possible that BMUs are simply initiated in damaged areas. Further, other stimuli clearly play a role, and thus a clear dichotomy between targeted vs. non-targeted views of BMU regulation is problematic and they need not necessarily be mutually exclusive (27). The classical view of remodeling has always been envisioned as a multi-functional role – including mechanical and physiological functions (26) such as calcium homeostasis. Efforts to disentangle the multi-faceted regulation of remodeling would be greatly aided by more and better 3D data regarding BMUs and related structures.

In sum, the capacity to directly test hypotheses related to regulation of BMU activity and/or the validation of *in silico* models is limited by a general lack of 3D data. Indeed, the activity of BMUs has largely been inferred from 2D observation of the secondary osteons they create. Our appreciation of the 3D structure of secondary osteons is similarly limited, and those data which are available (discussed below) consistently hint at greater structural complexity than commonly appreciated. Improving our 3D understanding of cortical bone microarchitecture would, thus, enhance our understanding of the remodeling process. As such, the objective of this paper is to provide an overview of methodologies for 3D investigation of cortical bone remodeling and, specifically, structures associated with BMU activity (resorption spaces) and the structures they create (secondary osteons). This review will survey a range of approaches spanning from histology to modern *ex vivo* imaging modalities, culminating with the growing potential of *in vivo* imaging. As such, it will span past, present, and emerging approaches. This collection of papers focuses on the theme of "putting the "*why*" back into bone architecture." Remodeling is one of two mechanisms "*how*" bone structure is dynamically modified, and thus an improved 3D understanding of this fundamental process is crucial to ultimately understanding the "*why*."

### Past: Histological Approaches

Long before the advent of modern 3D imaging modalities [e.g., confocal microscopy, scanning electron microscopy, and microcomputed tomography (micro-CT)], traditional light microscopy yielded a wealth of information about bone microarchitecture. Early applications of light microscopy revealed remodelingrelated structures including resorption spaces, mature osteons, and the canals within these osteons. The study of ground sections led to the first hypotheses pertaining to functional significance. Among the earliest of observations was a link between the extent of remodeling and age. Amprino and Bairati's (32) study of compact bone from humans and animals was one of the first to show an association between age and osteon population density (i.e. the number of osteons/mm2 ). Why this occurs, either through targeted replacement or through the accumulation of stochastic events, remains a central debate in the field. It is not surprising then than much remains to be learned from histology and that histology, in various forms, remains a mainstay of bone biology to this day. There are, however, limitations to extrapolating 3D structure from sections – a leap which requires idealization and assumption (33). Indeed, histomorphology can vary significantly due to differences in sampling location. The formation of an osteon by a BMU is an excellent example. Remodeling can create a new osteon anywhere along a bone's diaphysis and the resorption space in which it forms can extend over several millimeters (34). Stout et al. (35) found that sectioning osteons at different points along their lengths could result in the erroneous classification of different morphological "types." They provided the example of a "dumbbell-shaped" osteon being explained as the product of sectioning through a branching event. Beyond morphology of individual osteons, the rate of remodeling can vary within a bone. Skedros et al. (18), for example, demonstrated that the number of osteons present in the forelimb bones of Rocky Mountain mule deer increased in a proximal to distal fashion. While larger fields of view (FOV) have become increasingly feasible through the use of tiling microscopes, FOV has been a significant limitation of microscopy. This is especially problematic for the measure of rate of remodeling – reflected by the measure "activation frequency," which is defined as the number of new BMUs present in a particular unit of bone (10). Activation frequency can directly affect overall bone mass, because "remodeling space" temporarily reduces overall bone mass until it is subsequently filled in Ref. (36). Thus, a high remodeling rate can, in part, account for low mass. Activation frequency can be measured in either 2D sections (BMUs/mm2 /unit time) (10) or 3D volumes (BMUs/mm3 /unit time) based on the number of BMUs created per unit volume of cortex per unit time (37). The 3D extent of BMUs – their length or "range" – is a very difficult parameter to assess from histology, even from longitudinal sections. This parameter has the potential to impact the regional assessment of activation frequency with longer resorption spaces having increased potential to be detected in section. Cooper et al. (34) micro-CT analysis of 99 BMUs in human bone provided a measure of range varying between ~0.8 and 5.4 mm and most likely longer as they were limited to a 7 mm FOV and some remodeling events extended beyond. Keeping in mind the significant impact of sample site selection, activation frequency measures based on how *many* BMU-related resorption spaces are visible in sections acquired only millimeters apart from one another could be very different. Another problem with typical activation frequency measures is that, in section it can be difficult to differentiate BMU-related resorption spaces from other structures such as Volkmann's canals (36). That being said histology-based analytical techniques are still the gold standard for visualization of bone microstructure. Indeed, histology remains the primary means of investigating the relation between remodeling and microdamage. When identifying microdamage, in the form of microcracks, it is critical to ensure one is detecting *in vivo* damage and not artifacts of processing – particularly when utilizing grinding-based approaches as they can produce artifacts within bone that resemble cracks induced *in vivo* (38, 39). Microcracks are generally observed in 2D sections (40, 41), and thus their overall morphology, and relation to remodeling events, can be difficult to ascertain. For example, comparison of microcrack morphological features, such as size and shape, viewed in 2D has been shown to exhibit differences, to the extent that individual cracks look like entirely distinct structures (41, 42). This was seen to be the case in a study by Voide et al. (43), which found certain microcracks appeared linear and confined in the cross-sectional plane, whereas viewed in the perpendicular plane, they appeared diffuse.

Serial sectioning can alleviate some of these issues associated with 2D histology. Multiple sections increase the amount of bone analyzed and provide direct insight into the 3D nature of microarchitecture. Cohen and Harris' (44) analysis of serial decalcified sections from canine femora reported that osteons follow a spiral orientation around both the axis of the bone and axes of the osteons themselves. This seminal 3D study also reported that morphological characteristics vary in their distribution and this variation is linked to the specific locations within the bone. They observed that osteon cross-sectional area increased as they coursed distally throughout the bone and were also greater in the endosteal as opposed to subperisosteal regions (44), a finding suggestive of a mechanical influence on osteon morphology and thus BMU activity. Tappen's (45) work with block section staining using silver nitrate revealed morphological characteristics of osteonal canals, lacunae, osteocytes, and canaliculi. Tappen (45) highlighted numerous features of osteons, and the remodeling process which included levels of stain uptake in bone can be used as an indicator of mineralization; BMU-related resorption spaces are continuous with canals; these spaces can tunnel in opposite directions from one another or in only one direction, and some unforeseen force(s) dictate osteon/resorption space diameter as evidenced by the areas of mineralized bone resorbed in active areas of remodeling.

Despite its advantages, serial sectioning is challenging and thus only a few significant attempts to investigate cortical bone microstructure by this approach have been published (35, 44–46). Notably, these few studies have not found consensus with respect to some general aspects of osteon morphology. For example, Stout et al. (35) digital analysis of Tappen's (45) serial sections did not find evidence that osteons follow a spiral pattern around the axes of bone as suggested by Cohen and Harris (44). Rather, they found a complex pattern of branching osteons. Indeed, all 3D studies have revealed greater structural complexity than is frequently appreciated in the literature. While the literature is dominated by the "classical" BMU morphology presented in **Figure 1**, it should be noted that data indicating a complex diversity of remodeling-event types was reported in the mid 1960s by Johnson (4). We would argue that a lack of practical and efficient 3D techniques has contributed to the corresponding lack of information regarding BMU and osteon morphology. Beyond the difficult and tedious nature of existing serial-sectioning approaches, no small measure of luck is required to target these techniques to study individual BMUs and/or their related structures. Methodological difficulties ultimately led to the exploration of imaging-based approaches to the study of bone remodeling.

### Present: *Ex Vivo* Imaging

Since its introduction by Feldkamp et al. (47), X-ray micro-CT has become the "gold standard" for the non-destructive 3D analysis of trabecular bone morphology. More recently, this technology has increasingly been applied to cortical bone to analyze the porous network of canals in 3D (6, 14, 34, 48–53). With respect to the study of remodeling, micro-CT represents an excellent approach for detection of BMU-related resorption spaces since their distinctive cutting cones generally stand out from the relatively smaller canals of completed osteons. Another key advantage of micro-CT is its ability to non-destructively survey a large volume of bone for active remodeling events. An excellent example of this is our group's recent use of a laboratory micro-CT system (SkyScan 1172, Bruker, Belgium) to locate and trace BMUs in complete diaphyses of *Ursus americanus* (black bear) metacarpals and metatarsals without the requirement of sectioning (**Figure 2**) (54). In a more general sense, micro-CT can provide insight into remodeling through assessment of canal (and hence osteon) orientation (49) as well as the overall complexity of the canal network, which has been hypothesized to increase with cumulative remodeling activity (34, 48). The potential exists to examine osteon "steering" through assessment of localized canal orientation. Thus, the capabilities of micro-CT can alleviate or, at minimum, limit many issues associated with histological techniques for cortical bone: (1) interpolation of 3D structure from 2D sections is not needed as structures are observed directly and measured in 3D; (2) measures can be calculated in units of volume as opposed to area; (3) sample preparation is limited, eliminating structural alterations; (4) sample site location is less of an issue as visualization and analysis of larger volumes of interest – even entire bones – become feasible. These advantages are just beginning to be brought to bear on the *ex vivo* study of cortical remodeling, and micro-CT clearly has great potential to provide new insights through improved measures of morphology (including orientation, possible "steering," and 3D range) and activation frequency.

Despite its many benefits, micro-CT has several limitations. For cortical bone, the most notable is that while some delineation of osteon borders is possible from laboratory systems (55), this technology remains largely limited to detecting porosity at the vascular level and larger (e.g., osteonal canals and resorption

spaces). Synchrotron radiation (SR) micro-CT, with its many advantages conveyed by greater X-ray flux, its monochromatic coherence, and high brilliance as a result of a small source size can be utilized to visualize smaller-scale structures including osteon borders (56, 57) and different osteon morphologies (57), osteocyte lacunae (58–60) (**Figure 3**), and even microcracks (42, 43, 61, 62).

superimposed over BMUs. Diaphysis length **=** 31.84 mm.

A limitation shared by both micro-CT and SR micro-CT is that, in general terms, higher resolution comes at the cost of field of view. A further related limitation is that as resolution increases, so does the radiation dose. Improving resolution by a factor of 2 requires a dose increase by a factor of 16 to maintain image quality in micro-CT (63). This relation provides the practical limit on the resolution for *in vivo* X-ray based computed tomography. Thus, despite the fact that micro-CT imaging of trabecular bone microstructure in living animals is now commonplace (10–20 um voxel size), its application to the internal microstructure of cortical bone is not. An exception to this is the increasing use of high resolution-peripheral computed tomography (HR-pQCT) to analyze cortical porosity in humans, which will be discussed in the next section. The limited applications of micro-CT to the internal microstructure of cortical bone (e.g., osteonal canals) in animal models have all been *ex vivo* (64, 65). Working *ex vivo*, it is even possible to image microcracks utilizing barium sulfate as a contrast agent (66–68). With the vascular porosity of cortical bone beyond the reach of conventional *in vivo* micro-CT, it goes without saying that the resolutions required for microcrack imaging are not compatible with live animal imaging. Voide et al. (43) used SR micro-CT with a nominal

BMUs and osteocyte lacunae acquired by SR micro-CT at a 1.47 µm resolution (58). Image provided by Dr. Yasmin Carter.

resolution of 700 nm to visualize microcracks, *ex vivo*, in the femora of mice. Even if radiation dose was not an issue, holding a living animal still for sub-micron imaging would present a significant challenge in and of itself. Thus, while much can be learned about remodeling through *ex vivo* 3D imaging, *in vivo* detection and tracking of remodeling-related structures would represent a significant step forward with still greater potential to improve our understanding of this process. As will be discussed, there are specific advantages of synchrotron-based imaging which have created the potential for *in vivo* imaging of cortical porosity – opportunities that are just beginning to be explored and thus lie on the future horizon.

### Future: *In vivo* Imaging

A new generation of clinical research HR-pQCT scanners have given rise to a new opportunity for 3D, *in vivo*, characterization of human bone, both trabecular and cortical. With an isotropic voxel size of 82 μm, HR-pQCT has been at the forefront for microarchitectural analysis of the human appendicular skeleton (i.e., distal radius and tibia) (69, 71), as it is particularly advantageous for measuring cortical porosity (70). While HR-pQCT has been instrumental in studies concerned with osteoporosis induced change in microstructural cortical bone (71), its viability as a tool for targeting and tracking specific cortical structures, such as individual BMUs, is unclear due to resolution and other limitations associated with movement artifacts and the restriction to primary trabecular bone sites (i.e., wrists and ankles) (72). When one considers that the size of a single HR-pQCT voxel is on the same scale as the average canal diameter within human cortical bone, ~120 μm for pooled sexes (6), it becomes clear that this technology cannot resolve all osteonal canals. A recent study comparing HR-pQCT as a tool for measuring pore sizes against SR micro-CT showed that the accuracy of HR-pQCT was a factor of the size of the pores measured (70) – current resolution limits its ability to recognize small pores, thus introducing partial volume effects into analyses (72). That said, as BMU-related resorption spaces average 250 μm in diameter [e.g., the average diameter of an osteon in humans (22)], HR-pQCT may indeed hold the potential of tracking individual remodeling events in human cortical bone. Thus, while HR-pQCT represents a powerful new tool for the assessment of cortical porosity, it is limited to the largest of cortical pores located in the ultra-distal peripheral skeleton.

Human-based studies are clearly an important avenue for study, but more mechanistic studies will require controlled model systems. When considering animal models – particularly smaller animals – the punishing relation between resolution and radiation dose complicates *in vivo* imaging. The risks associated with high X-ray doses are not only harmful to a living subject in an acute sense, they include the potential to alter bone structure, leading to osteopenia, growth arrest, fracture, and malignancy (73) through altered osteoclast and osteoblast activity (74). For example, at a voxel size of 10.5 μm and a radiation dose of 845.9 mGy, longitudinal scans of the hind limb in mice acquired over a 5-week period resulted in decreased bone volume and increased trabecular separation (75). Even though SR micro-CT makes *in vivo* imaging of cortical porosity more plausible through higher resolution and shorter scan times, the limitations imposed by radiation dose are still a factor. As will be discussed, the monochromatic (single X-ray energy) capabilities of synchrotrons enable alternative methods of developing contrast beyond absorption. These so-called "phase contrast" techniques are opening the door to new possibilities for *in vivo* imaging.

Several forms of phase contrast imaging, including in-line phase contrast, diffraction-enhanced imaging, and interferometer-based imaging, have been increasingly explored for their potential utilization in biomedical applications. An excellent overview of these methodologies is provided by Zhou and Brahme (76). Unlike conventional radiography, where images are based on the differences in X-ray absorption related to an object's internal structure, phase contrast images contain information related to the refractive index of an object because the images produced are derived from the pattern of interference created from diffracted and undiffracted waves (77). Differences in refractive index of the sample's internal structure will refract or bend the X-ray wavefront as they pass through the target and these differences can be used to generate contrast (76). For in-line phase contrast imaging – the simplest of these techniques to implement – detection of changes in X-ray refractive indices within a sample can be optimized by altering the distance between the object and the detector (76, 78) (**Figure 4**). Simple implementation and enhanced detection of internal features have made the fusion of SR computed tomography and phase contrast imaging particularly advantageous for visualization of microscopic bone structures such as osteocyte lacunae (58, 79, 80), and even nanoscopic structures such as the lacuno-canalicular network (59). Phase contrast is not dependent on X-ray attenuation and thus higher energies, where attenuation is lower, can be employed (76). This creates the possibility of higher resolution with equal or reduced dose compared with attenuation-based imaging.

Our group recently tested the feasibility of in-line phase contrast SR micro-CT for imaging the cortical porosity in the forelimb of rats with radiation doses comparable to those commonly employed *in vivo* for imaging trabecular microarchitecture (78). Since its' first implementation (81), numerous *in vivo* SR micro-CT protocols varying in both dose and image resolution have been utilized to visualize the microarchitecture of bone in small animal models (78, 82–84). Currently, there is no consensus for what is considered to be a safe dose for live animal *in vivo* imaging. Voxel sizes in the 10–15 μm range are typical for laboratory micro-CT with dose rates of 0.4, 0.338, and 0.939 Gy, respectively (85–87). However, at a voxel size of 11.7 μm, Matsumoto et al.'s (83) *in vivo* scans of the knee joint in living mice involving 5 Gy showed no signs of radiation sickness. Based upon this context, our goal was to create an *in vivo* imaging protocol involving a dose as close as possible to those used routinely in conventional micro-CT and no higher than 5 Gy.

FIGURE 4 | (A) Schematic of attenuation based X-ray imaging where images are produced based on the degree of absorption relative to an object's internal structure. (B) Schematic of in-line phase-contrast imaging based on an object's refractive properties. As X-rays target an object at different angles, variations in the object's internal structure will refract the X-rays and cause a shift in the light wave as it propagates through and increasing the object-todetector distance will produce a contrast image (76). Reprint permission granted by the publisher.

In a proof-of-principle study, *ex vivo* data collected at the biomedical imaging and therapy (BMIT) facility of the Canadian light source (CLS) synchrotron were compared against laboratory micro-CT protocols and their related doses. Studying rat forelimbs, it was found that SR provided superior detection of cortical pores without a substantial increase in dose (11.8 μm voxels, 2.53 Gy) beyond that used in the laboratory systems (18 μm voxels, 1.2–1.5 Gy; 9 μm voxels, 11.7–18.2 Gy) (78) (**Figure 5**). For in-line phase SR micro-CT, the optimal target-to-detector distance for our configuration was found to be 0.9 m. Subsequent to the *ex vivo* trials, the first *in vivo* trial was performed (78) and the first longitudinal studies employing this protocol are under way – providing encouraging results with respect to tracking individual remodeling events (**Figure 6**).

Synchrotron radiation micro-CT, taking advantage of in-line phase contrast, has brought about a new opportunity for *in vivo* imaging to directly test causative hypotheses relating to cortical bone remodeling; however, it too is not without its limitations. Access to synchrotrons is inherently limited due to the relative

FIGURE 5 | Reconstructed slices of rat forelimbs depicting visualization of cortical porosity based on the imaging system used: (A) *in vivo* laboratory SkyScan 1176 micro-CT (18 μm, 1.2–1.5 Gy dose), (B) *in vivo* synchrotron micro-CT slice measured using the C4742-56-12HR camera (11.8 μm, 2.53 Gy dose), (C) *in vivo* laboratory SkyScan 1176 (9 μm, 11.7–18.2 Gy dose) (78). Reprint permission granted by the publisher.

scarcity of such facilities and, in particular, those with biomedicalfocused facilities capable of imaging live animals. Even with the high resolution afforded by SR micro-CT and the ability to minimize dose through use of in-line phase contrast imaging, some structures such as osteon borders and microcracks are not observable, *in vivo.* Indeed, such structures remain challenging to image *ex vivo*, requiring very high-resolution systems. Movement artifacts are problematic as the difficulty in holding the animal perfectly still increases along with resolution. The smallest of the standard animal models (e.g., mice and rats) have the further drawback of not exhibiting much, if any, natural cortical remodeling – although it can be induced by experimental means [e.g., fatigue loading in the rat (31)]. Notably, this imaging approach has the potential to be extended to larger animal models including those such as the rabbit, which exhibit remodeling under normal conditions. Larger cortical canals in such animals could also relax the resolution needs and, concomitantly, reduce radiation dose.

## Conclusion

The remodeling process, carried out by the activity of BMUs, is of great interest to the bone biology community. This process represents the primary means of skeletal change after maturity and lies at the root of many chronic bone diseases, including osteoporosis. Visualization of cortical bone microarchitecture and, specifically, the remodeling process have progressed from 2D histological analysis through *ex vivo* 3D imaging and now to *in vivo* 3D analysis. This progression has paralleled but lagged behind visualization of trabecular bone microarchitecture due to the smaller scale of the target features and the need for high resolution which suffers the complication of increased radiation dose for X-ray based imaging. An important consideration is the caveat that while this progression involves a decreasing level of sample destruction/invasiveness, there is a trade off in terms of the structures one can visualize. The imaging approaches are best suited to detection of porosity. Thus, for some applications, histology (serial or otherwise) remains the most powerful or possibly only approach available. That said, there is a great potential to combine the strengths of these approaches – fusing 2D and 3D imaging to maximize the information available. This will enable targeted histology – directed by imaging data. Looking to the near future, we believe that this approach will see application in direct testing of hypotheses related to the regulation of BMU activity – including questions related to steering/orientation, the role of microcracks, and the relation to other stimuli, including possibly cellular signals. Such data will, in turn, prove invaluable for validating *in silico* models – an area of increasing focus in bone biology. Ultimately, we believe the novel insights possible through 3D data will shed significant new light on the "how" of bone aging, adaptation, and disease. Understanding the "how" is critical to understand the "why."

## Acknowledgments

The authors thank the organizers of this topic, Drs. Phil Salmon, Daniel Chappard, and Andrew Pitsillides. The authors would also like to thank: the Canadian Light Source BMIT line staff, particularly Drs. George Belev, Ning Zhu, and Adam Webb, Dr. Melanie van Der Loop, and staff of the University of Saskatchewan's Laboratory Animal Services Unit, and the support of Dr. Mike Doschak of the Pharmaceutical Orthopedic Research Lab at the University of Alberta, and the contributions of Dr. Dean Chapman, Dr. Yasmin Carter, Isaac Pratt, and Danielle Kabatoff. Finally, the authors thank Ms. Suyoko Tsukamoto for providing black bear bone samples (MB Conservation WBO8992). Support for this research was provided by the Natural Sciences and Engineering Research Council (NSERC) of Canada via a Discovery Grant to DMLC (RGPIN-2014-05563) and the Canadian Foundation for Innovation and Canada Research Chairs program. KH is supported by the University of Saskatchewan, College of Medicine, and is a CIHR-THRUST Fellow.

### References


**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Harrison and Cooper. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## 3D porous architecture of stacks of **β**-TcP granules compared with that of trabecular bone: a microcT, vector analysis, and compression study

*Daniel Chappard1,2\*, Lisa Terranova1 , Romain Mallet1,2 and Philippe Mercier1,3*

*1GEROM Groupe Etudes Remodelage Osseux et bioMatériaux – LHEA, IRIS-IBS Institut de Biologie en Santé, CHU d'Angers, L'Université Nantes Angers Le Mans, Angers, France, 2Service Commun d'Imagerie et Analyses Microscopiques (SCIAM), IRIS-IBS Institut de Biologie en Santé, CHU d'Angers, L'Université Nantes Angers Le Mans, Angers, France, <sup>3</sup> Laboratoire d'Anatomie, Faculté de Médecine, Angers, France*

### *Edited by:*

*Phil Salmon, Bruker microCT, Belgium*

#### *Reviewed by:*

*Egon Perilli, Flinders University, Australia Peter James Delisser, University of Bristol, UK*

*\*Correspondence: Daniel Chappard daniel.chappard@univ-angers.fr*

#### *Specialty section:*

*This article was submitted to Bone Research, a section of the journal Frontiers in Endocrinology*

*Received: 17 December 2014 Accepted: 28 September 2015 Published: 12 October 2015*

#### *Citation:*

*Chappard D, Terranova L, Mallet R and Mercier P (2015) 3D porous architecture of stacks of* β*-TCP granules compared with that of trabecular bone: a microCT, vector analysis, and compression study. Front. Endocrinol. 6:161. doi: 10.3389/fendo.2015.00161*

The 3D arrangement of porous granular biomaterials usable to fill bone defects has received little study. Granular biomaterials occupy 3D space when packed together in a manner that creates a porosity suitable for the invasion of vascular and bone cells. Granules of beta-tricalcium phosphate (β-TCP) were prepared with either 12.5 or 25 g of β-TCP powder in the same volume of slurry. When the granules were placed in a test tube, this produced 3D stacks with a high (HP) or low porosity (LP), respectively. Stacks of granules mimic the filling of a bone defect by a surgeon. The aim of this study was to compare the porosity of stacks of β-TCP granules with that of cores of trabecular bone. Biomechanical compression tests were done on the granules stacks. Bone cylinders were prepared from calf tibia plateau, constituted high-density (HD) blocks. Low-density (LD) blocks were harvested from aged cadaver tibias. Microcomputed tomography was used on the β-TCP granule stacks and the trabecular bone cores to determine porosity and specific surface. A vector-projection algorithm was used to image porosity employing a frontal plane image, which was constructed line by line from all images of a microCT stack. Stacks of HP granules had porosity (75.3 ± 0.4%) and fractal lacunarity (0.043 ± 0.007) intermediate between that of HD (respectively 69.1 ± 6.4%, *p* < 0.05 and 0.087 ± 0.045, *p* < 0.05) and LD bones (respectively 88.8 ± 1.57% and 0.037 ± 0.014), but exhibited a higher surface density (5.56 ± 0.11 mm2 /mm3 vs. 2.06 ± 0.26 for LD, *p* < 0.05). LP granular arrangements created large pores coexisting with dense areas of material. Frontal plane analysis evidenced a more regular arrangement of β-TCP granules than bone trabecule. Stacks of HP granules represent a scaffold that resembles trabecular bone in its porous microarchitecture.

Keywords: porosity, microCT, 3D packing, 3D geometry, bone graft, fractal, β-TCP, granules

## INTRODUCTION

Bone is a calcified connective tissue, which fulfils a large number of functions in the body (1). Being both rigid and elastic due to the unique composition of its matrix made of collagen type I and hydroxyapatite crystals, it is adapted to gravity and muscle contractions. Bone cells are able to remodel and adapt bone mass and architecture throughout life allowing repair of microdamage (microcracks, trabecular microfractures) or fracture. To reduce the amount of this material within the body, Nature has developed a number of strategies to provide the maximum biomechanical competence with a genetically controlled bone mass. During aging, epigenetic factors, such as gravity and muscle strains, also modify bone mass and microarchitecture (Wolff 's Law). Trabecular microarchitecture is now a well-recognized type of self-organization that fulfils this condition (2, 3).

Bone loss can occur systemically, for example, during aging. However, localized bone loss can also occur, as in some intra-bony defects, such as solitary or aneurysmal cysts, or acquired bone defects as observed after tooth extraction or during hip prosthesis revision. Granular biomaterials are often used to fill these defects, especially in non-bearing bones. Among the different orthophosphates that can be produced in large quantity by industry, beta-tricalcium phosphate (β-TCP) appears the most suitable, particularly in maxillo-facial and dental surgery (4). Granules of 1000–2000 μm diameter represent the most commonly used size for filling alveolar sockets or increasing bone volume by sinus lift (5). However, when the surgeons place granules within a bone defect, it is the voids between the granules that represent the interconnected space available for vascular sprouts and bone cells to invade the grafted area. Calcium phosphate ceramics are known to be too brittle to be placed in load-bearing areas. Ceramic granules allow a perfect filling of the bone defects, but they are limited to grafted areas where no significant strength is needed (alveolar ridge augmentation, sinus lift …) (5). Little is known about the 3D arrangement of a grafted stack of granules, and if the porosity of the stack is similar to that of trabecular bone to allow the accessibility of biological fluids, body fluids, cells, and vascularization within the grafted site (6, 7). In addition, when filling a bone defect, the surgeon has to avoid crushing of the granules to avoid altering the granules' microarchitecture (8). Hence, the biomechanical resistance of the granules in compression is also an important factor to be taken into account when choosing the grafting material (9).

The aims of the present study were (1) to compare the pore microarchitecture of the two types of β-TCP granules providing stacks of low porosity (LP) and high porosity (HP), respectively, with that of the two types of trabecular bone [high density (HD) and low density (LD)], in order to assess the granular stack that had the pore microarchitecture more similar to trabecular bone; (2) to compare the biomechanical properties in compression of these two types of granules."

### MATERIALS AND METHODS

### Bone Samples

For the bone with HD, five bone cylinders were harvested from the upper tibia epiphysis of young butchered animals (calves, aged around 15 months, 500–600 kg). The calibrated cylinders (10 mm in diameter) were obtained with a diamond-tipped coring tool (Starlite Industries, Rosemont, PA, USA) perpendicular to the tibial plateau. Similarly, five cylinders (10 mm diameter) with low bone density (LD) were prepared from the upper tibial epiphysis from five aged human cadavers from the Anatomical Laboratory (mean age 78.3 ± 5.2 years). For both the HD and LD cylinders, the articular and subchondral bone were removed to provide cylinders 23–25 mm in length containing only the secondary spongiosa. The subjects had given their body to science before death and had completed a form indicating that their corpses were to be used for medical education and research. Subjects were examined upon arrival at the laboratory of anatomy of the Faculty of Medicine. All bone samples were harvested following a strict protocol of the laboratory of anatomy; e.g., it was impossible to know any detail which could lead to the subject's identification (causes of death, previous medical records).

## **β**-TCP Granules

Beta-tricalcium phosphate granules were obtained from Kasios (L'Union, France). Granules were prepared by polyurethane foam technology and used as received. This technique has been extensively described elsewhere (10, 11). Two types of granules were used and placed in polyethylene test tubes (10 mm in diameter): compact granules were prepared with 25 g of β-TCP grain in a fixed volume of slurry and provided a LP stacks; lighter granules with 12.5 g of β-TCP grain 1 in the same volume of slurry produced the HP stacks (12).

The morphology of the two types of granules was qualitatively evaluated by scanning electron microscopy (SEM, Jeol JSM 6301F) before analysis. The β-TCP granules were glued on an aluminum stub with a double-sided carbon sellotape. They were gold-coated by sputtering with an MED 020 (Bal-Tec, Chatillon sur Cher, France). Images were captured at a 3-kV acceleration voltage in the secondary electron mode.

### Biomechanical Analysis

The compression characteristics of the two types of β-TCP granules were determined on an Instron 5942 machine with Bluehill software (Instron, Elancourt, France) (**Figure 1A**). As the test is destructive, new intact granules were placed into the cylindrical test cell of the machine and compacted with an adapted punch 10 mm in diameter. The maximum applied load was fixed at 400 N with a 2-mm min<sup>−</sup><sup>1</sup> displacement. For each type of stacks of granule (HP or LP), the load–displacement curve was recorded (N/mm) automatically by the machine. Stiffness was automatically provided by measuring the slope of the tangent to the curve at the beginning of the deformation (N/mm). The maximum displacement of the punch (in mm) was obtained when the load reached 400 N. The work to failure (N⋅mm) was determined as the area under the load–displacement curve from yield until failure and represents the energy absorbed by bone during plastic deformation (**Figure 1B**) (13). These parameters were calculated according to previously published equations (14, 15). Analysis was repeated five times for each type of granules and the averaged values were considered.

## MicroCT Analysis

MicroCT analysis was performed on a Skyscan 1172 microcomputed tomograph (Bruker microCT, Kontich, Belgium). Bone samples were placed in polystyrene test tubes (10 mm in diameter) and scanned while in fixative (10% formalin). Granules of β-TCP were placed in similar tubes and gently agitated to allow the granules to settle and optimize their 3D spatial arrangement in the dry state. Tubes were filled with β-TCP granules up to 30 mm in height. For each series of granules, analysis was done in triplicate. The microCT was operated at 80 kV, 100 μA with a 0.5-mm aluminum filter. The pixel size was fixed at 4.95 μm and the angular rotation step at 0.25°. 3D models of bones and granules were obtained with VG Studiomax (Volume Graphics GmbH, Heidelberg, Germany) operating in the volume rendering mode. Microarchitecture parameters were obtained with the CTAn software (Bruker) after global thresholding of the pores (16). They included measures of porosity (Po, in %), mean pore diameter (Po⋅Diam, in μm), and specific surface [representing the interface between bone or biomaterial and the porosity, BS/

TV, in mm2 /mm3 , where BS is the whole pore surface and TV is the volume of interest according to the ASBMR nomenclature (17)]. For the bone cylinders and the granule stacks, the region of interest (ROI) over-imposed on the 2D binarized sections was a 5-mm wide square; the height of analysis of the samples has concerned 15 mm. Pore diameter was determined by the spherefitting method in the CTAn software (18). In order to further investigate the porosity, the binarized stacks of 2D microCT sections of each sample were analyzed by software written in Matlab (MathWorks, Natick, MA, USA). The algorithm for evaluation of porosity used a vector projection on a frontal plane (Vectopor) and has been extensively described elsewhere (19). In the present study, the pixels of the ROI which belonged to the same column received the same pseudo-color according to the number of pixels superimposed on pores. A look-up table (LUT) was designed ranging from deep blue (no porosity) to red (high porosity). Intermediate values depend on the value of porosity along the vector according to the LUT. The frontal plane vector projection image was assigned intensity values for each pixel based on the number of pixels, which contained pores in the related column of the microCT binarized image. The frontal plane image was saved with the colorized LUT and analyzed by the FracLac plugin developed for ImageJ (20, 21). The box-plot fractal dimension (*D*f) and the lacunarity (λ) of the frontal plane images were determined as previously described (19). The boxplot fractal dimension evaluates the complexity of the material to fill the reference space and is helpful in the study of porous objects, such as bone (22, 23) or materials (24, 25). Lacunarity is another recently described fractal parameter reported to improve the description of a fractal porous object (26, 27). Lacunarity is influenced by the variation of the pores in a given structure; a low lacunarity reflects homogeneity, while a high lacunarity is an indicator of heterogeneity (28, 29).

### Statistical Analysis

Statistical analysis was done with Systat 13 (Systat Inc., San José, CA, USA). For each parameter, results are expressed as mean ± SD. Differences for porosity between the four groups (HD, LD, HP, and LP) were sought with the Kruskal–Wallis nonparametric analysis of variance followed by the Conover–Inman test for all pairwise comparisons. For biomechanical parameters, differences between the HP vs. LP groups were analyzed by the non-parametric Mann–Whitney *U* test. A difference was considered as significant when *p* < 0.05.

### RESULTS

### Qualitative Examination of the **β**-TCP Granules

The difference in morphology of the granules was evident to the naked eye and confirmed by SEM. The 12.5-g preparation of granules was the most porous with thin walls of sintered β-TCP, while the 25-g granule preparation was denser, with numerous concave surfaces and occasional macropores (**Figure 2**). In all cases, the inner porosity due to sublimation of the polyurethane foam was clearly identified as triangular voids connected by thin channels.

### Biomechanical Analysis of the **β**-TCP Granules

The load vs. displacement curves of both HP and LP types of granules are shown in **Figure 1C**. The first domain of the curve corresponds to the collapse of the granules, after which the deformation is due to their irreversible compression. The HP stack had a statistically significant higher maximum

granules prepared with 25 g in the slurry, (B) HP granules prepared with 12.5 g. Note the inner porosity due to the sublimation of the polyurethane foam as triangular voids connected by thin channels (arrows).

displacement compared to the LP stack, as expected (*p* < 0.05, **Table 1**). On the contrary, the stiffness and the work to failure were significantly higher in the LP stacks compared to the HP stack (*p* < 0.05 in both cases).

### Morphometric Analysis of Bone and Granules

A 3D-rendering of the two types of bone cores and granule stacks is shown in **Figure 3**. Morphometric parameters are listed in **Table 1** and illustrated in **Figure 4**. The HD bone cores had a significantly smaller porosity than LD bones. The LP stacks had the lowest porosity among all groups. Conversely, HP stacks presented a porosity that is intermediate between that of a HD bone and that of a low-density bone. The HP formulation produced a ~28% increase of porosity associated with a reduction of the mass of β-TCP. This was associated with a significant increase in BS/TV in the HP group compared to the LP group (+11.6%). The vector analysis of porosity clearly evidenced differences between the two pairs of materials: HD bones appeared inhomogeneous with areas containing large amounts of bony material (blue areas according to the LUT), while LD had a large number of red spots indicating HP (**Figures 5A,B**). Very similar findings were observed for the LP and HP stacks of granules: LP stacks had large areas of material but large holes were also present. On the contrary, the morphology of the stacks of HP granules had a more regular distribution of the pores and material (**Figures 5C,D**). Interestingly, the fractal dimension *D*f did not differ among the four groups, but lacunarity λ was high in HD bones and even higher in the LP granules (**Table 1**).

### DISCUSSION

In this series of materials comprising bone cores and granular stacks made of β-TCP granules, a quantitative and morphological analysis of porosity was done. Porosity of the LP stacks was noticeably lower than a HD trabecular bone. However, the porosity of the HP stacks (and conversely the amount of material) was in the range of normal bone volume as determined by classical


TABLE 1 | Morphometric parameters of bone cylinders with a high-density (HD), low-density (LD), and **β**-TCP granules (compact low porosity granules: LP; high porosity granules: HP).

*Data are provided as mean* ± *SD. Differences between the four groups were sought with the Kruskal–Wallis analysis of variance followed by the Conover–Inman post hoc test for all pairwise comparisons. For biomechanical parameters, differences between the HP vs. LP groups were analyzed by the Mann–Whitney U test.*

*\*Significantly different from HD cylinders with p* < *0.05. § Significantly different from LP granules with p* < *0.05.*

Frontiers in Endocrinology | www.frontiersin.org 78 October 2015 | Volume 6 | Article 161

FIGURE 3 | Microtomographic analysis of the different types of materials (appearing in white). (A) High-density bone (HD) with numerous trabeculae and low porosity; (B) low-density bone (LD) from aged subjects with reduced trabeculae and increased porosity; (C) of compact granules (LP) of β-TCP; (D) scaffold of HP granules. All images have been prepared with a square ROI of 5 mm side; so the volume of interest is a parallelepiped (height: 15 mm; square section: 5 mm).

histomorphometry methods (30). This reflects the microarchitecture of the stacks of HP granules which contain more macropores and thin branches than the LP granules. However, as expected, this effect is detrimental to the biomechanical properties of the HP granules whose stacks have lower compressive mechanical properties, as evidenced by lower stiffness values, lower work to failure, and maximum displacement. The surgeon should be aware not to crush the granules in the grafted areas and the use of a bone syringe with a large diameter may help.

It should be noted that the BS/TV of HP granules is significantly higher than that of LP granules or that of HD and LD bones. This is especially important since BS/BV represents the bioactive surfaces of β-TCP that will be in direct contact with the biological fluids and the cells (31). The BS/TV values obtained in this study with β-TCP granules are similar to those reported by other groups with BCP (a mixture of β-TCP and hydroxyapatite) (32). Recent observations have stressed the interest of increasing the specific surface (BS/TV) of orthophosphate materials to promote

bone formation (33). So, the HP formulation represents the highest surface available for osteoblast attachment and proliferation (osteoblasts are bone-forming cells). Because β-TCP is a biomaterial colonized by osteoconduction from nearby trabeculae, this parameter is of great importance (34).

MicroCT analysis is a relatively recently introduced tool in the analysis of bone and materials. Values for parameters, such as BV/TV and Po, are highly correlated with histomorphometry on 2D sections (35). In the present study, the difference in porosity between the stacks of granules can be clearly identified as previously reported (12). Different types of parameters have been proposed to evaluate the 3D arrangement of porous materials (e.g., the box-plot fractal dimension) (22, 23) and the interconnectivity of the pores created between the granules (29). The method was successfully applied to evaluate the complexity of different types of granule aggregates: sandstone (36), coal (37), metal grains (38), and soils (39). The complexity of the granule repartition, evaluated by *D*f, did not differ between the four groups. However, fractal lacunarity appears as an interesting new tool capable to differentiate the groups but it does not bring morphological information (12, 28, 29, 40). We have recently developed vector analysis for visualization and quantification of the porosity in objects presenting a complex form (19). The method was found useful to analyze pores in materials prepared with various types

minimal porosity in blue according to the LUT (at the bottom of the image). All images have been prepared with a square ROI of 5 mm side length and the scale bar represents 1 mm.

of porogens and in bones having a particularly complex shape. This is the case with alveolar bone at the mandible, which is very difficult to analyze because of the presence of the tooth roots. Recently, the bone loss induced by paralysis of *Mus masseter* and *Mus temporalis* was evidenced by this technique (41). The frontal plane is produced line by line during the vector analysis of each microCT image and represents the projection of porosity. The complexity of the image can be measured by fractal geometry on this projected 2D image. In the present series, there was no difference for *D*f among the four groups under study. However, lacunarity is known to be influenced by the uniformity of the

### REFERENCES


spatial distribution; it is a useful parameter when objects do not differ by their fractal dimension (42, 43). The maximal λ value was reached in the LP (25 g) group meaning that very large pores are created in the 3D arrangement of the granules. On the other hand, the lowest λ values were obtained for the HP (12.5 g group) granules and the LD bone, which presented a regular distribution of pore dimensions. Similar relationships have been reported in studies concerned with the porosity of soils (44, 45).

The presence of pores larger than 100 μm has been recognized as a key factor in the vascular invasion of grafted materials into bone (46, 47). In this study, Po⋅Diam values for bone were within the values typically reported for bone (48); Po⋅Diam of both HP and LP are larger than 300 μm and well suited for body fluid and cell invasion (49).

There are some limitations to this study: (1) no compressive analysis was done on the bone cylinders, but synthetic ceramics are well known to be more brittle than bone; (2) bovine bone was used to obtain the HD cylinders, but the biomechanical parameters for trabecular bone are very close (stiffness in compression for bovine bone: 173 MPa vs. human: 139–472 MPa) (50); and (3) the other microarchitectural characteristics of the β-TCP materials having been reported elsewhere are not duplicated in the present study (12).

## CONCLUSION

The present study shows that the two types of β-TCP granules prepared differently did not produce the same interconnected porosity: LP granules provided very large but less numerous pores and HP granules provided a HP with a large interface. The HP stacks of granules presented a porosity similar to trabecular bone, although the granules were physically independent. This study confirms that vector analysis is a suitable method for the evaluation of porosity of complex objects.

### ACKNOWLEDGMENTS

The authors thank Mrs. Laurence Lechat for secretarial assistance. This work was made possible by grants from ANR, program LabCom "NextBone." The authors thank Kasios, 18 chemin de la Violette, 31240 L'Union – France for providing the different formulations of granules. The analysis software (VECTOPOR) is now licensed by the Agence de Protection des Programmes (APP). Authors also thank Mrs. Nadine Gaborit for her skillful technical assistance with microCT and Phil Salmon for having reviewed the manuscript.

in aging and in osteoporosis. Implications for the microanatomic and cellular mechanisms of bone loss. *J Clin Invest* (1983) **72**:1396–409. doi:10.1172/ JCI111096


macroporosity, percentage, on, bone, ingrowth. *Biomaterials* (1998) **19**:133–9. doi:10.1016/S0142-9612(97)00180-4


**Conflict of Interest Statement:** The authors declare that this research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2015 Chappard, Terranova, Mallet and Mercier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*