Cristian Enachescu
email: cristian.enachescu@uaic.ro

Roman Stoleriu

Masamichi Nishino

Seiji Miyashita

Roman Stancu

Maciej Lorenc

Roman Bertoni

Hervé Cailleau

Eric Collet

Laurentiu Stoleriu

Alexandru Stancu

A theoretical approach for elastically driven

Keywords: numbers: 75.30.Wx, 64.60.-i, 75.60.-d

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Spin crossover inorganic compounds are among the few systems showing bistability at a molecular scale with switchable properties by variations of temperature and pressure or by light irradiation [START_REF] Gütlich | Spin crossover in transition metal compounds[END_REF][START_REF] Halcrow | Spin-crossover materials -properties and applications[END_REF] . Pulsed laser irradiation opens the possibility to control this bistability in an ultrafast way, through strongly out-of-equilibrium conditions [START_REF] Collet | [END_REF][4][5][6][7][8] . Their molecular switching occurs between the paramagnetic low spin state (LS), stable at low temperatures and paramagnetic high spin (HS) state of higher entropy, stable at higher temperature. Remarkably, not only magnetic and optical properties are different in the two states, but also the molecular volume, being around 3% higher in the HS state due the occupancy of antibonding ligand orbitals in HS molecules 9 . This size changes during thermal conversion between LS and HS states induces elastic interactions mediated by the lattice, which, if high enough, are responsible for abrupt transitions and cooperativity phenomena. At low temperatures, such crystals can be photo-switched from the ground LS state to the long-lifetime metastable HS state. By exciting some specific charge transfer bands, with a broad band lamp or a laser, the so called Light Induced Excited Spin State Trapping -LIESST effectoccurs 10,11 . At long time scales, varying from nanosecond to days as a function of compound and temperatures, the metastable HS state relaxes towards the ground state by a nonadiabatic, nonradiative process, either via tunnelling -at low temperature -or thermally activated process 12 . The recent experiments using femtosecond laser for photoexcitation and the detailed analysis of the subsequent evolution of the system allowed the study of very fast dynamical phenomena and contributed to a better understanding of the behaviour of spin crossover compounds at molecular level 5,13,14 . Observations of molecules in solutions or in solid compounds demonstrated that the photoswitching from the LS towards the HS states, involving both change of electronic states and molecular structural reorganization occurs within ≈ 200 fs [15][16][17][18] . In solids, this initial photo-switching of few molecules from LS to HS states triggers further increases of the fraction of HS molecules X HS . This elastically driven cooperative switching of spin crossover compounds occurs due to the propagation of the elastic distortions in the medium, the superposition of which induces volume expansion of the lattice. This so-called elastic step 5 during the of the out-of-equilibrium response following fs laser excitation, takes place on the acoustic time scale, which corresponds to the ratio between relevant size scale of the sample and the sound velocity in the medium. The photoresponse measured through the time evolution of the fraction X HS of HS molecules, is directly probed by ultrafast spectroscopy. Figure 1 shows typical photoresponse to fs laser excitation of SCO nanocrystals of

[ F e III (3 -M eO -SalEen) ] P F 6 , which was reported in detail in ref. 5 . It is non-linear with excitation density. The fraction of molecules switched by light from LS to HS states X hν HS , measured at the picosecond time scale, depends linearly on the excitation density because of the local nature of the photoswitching process and each photon switches a single molecule. When the fraction of photoswitched molecules from LS to HS state is below 3%, the LS ground state is recovered within ∼ 1 µs.

At higher excitation density the HS fraction increases. For high excitation density, the molecular swelling overcome the lattice pressure and drives its expansion, favouring so the HS state of higher volume. In this case a self-amplified transformation is observed and, in average, up to 7 molecules are switched per photon. The additional conversion on the elastic step is then defined as ∆X HS = X El HS -X hν HS . Therefore below the amplification threshold ∆X HS < 0, whereas above ∆X HS > 0.

Among the many models elaborated for the study of spin crossover compounds, elastic models [19][20][21][22][23] are the more suitable for the study of the elastic step, as they build on the common sense by which interactions between molecules are of elastic origin and stem from the size difference between HS and LS molecules. Here we use a dynamic mechanoelastic model [23][24][25] applied for open boundary systems in order to emulate the evolution of the system after femtosecond excitation as function of several parameters. Ab-initio models are useful for study the details of electronic and bonding states 26,27 , but for studies of phase change, the mechano-elastic model and other spin-crossover elastic models have advantage. Especially for the phenomena of elastic step, the combination of molecular dynamics and stochastic dynamics, treated in this study, is effective. This model was briefly introduced in ref. 5 and in the present paper we discuss the role of different parameters in the model and the resulting dynamics due to two competing and opposite going processes, namely the local relaxation and the non-local amplification.

II. THE MODEL

In the mechanoelastic model, the spin-crossover molecules are represented as rigid spheres, situated in triangular configurations in open boundary lattices and interacting through connecting springs, with the elastic constant k. The conversion between the two states, HS and LS, is accompanied by a variation of the volume of a switching molecule, producing an instantaneous elastic force in its neighboring springs, which is transferred, step by step, to all the molecules in the system, in a degree depending on their relative position form the switched molecule. Therefore, the volume and the shape of the whole system dynamically change with every molecular switch, as elastic deformation propagates through the sample.

The motion of molecules is determined by local elastic forces and is accounted by solving the following differential equations for all molecules in the system, for a specific number of steps n 1 .

       m d 2 x i dt 2 = F i,x -µ dx i dt m d 2 y i dt 2 = F i,y -µ dy i dt , ( 1
)
where x i , y i are Cartesian coordinates of the molecule i, m is the mass of the molecule, µ is a damping constant, and F i,x , F i,y the components of the instantaneous force acting on molecule i given by the sum of the forces from the neighbouring springs:

           F i,x = ∑ neighbor springs kδr ij,x F i,y = ∑ neighbor springs kδr ij,y , (2)
where δr ij is the deviation of the neighbour spring j from its unstressed value. In order to study dynamical properties of the system during a spin transition within the ME model, one assigns for every spincrossover molecule a HS → LS and a LS → HS transition probability, depending on the temperature T , on intrinsic molecular parameters (the HS -LS energy difference D, the degeneracy ratio g, the effective activation energy E A ) and on the elastic interactions between molecules, represented here by the local pressure p i acting on molecule i, defined as

p i = ∑ neighbor springs k A δr ij , with A
the molecular cross-sectional area. These Arrhenius-type probabilities can be written as:

P i HS→LS = 1 τ exp ( - E A -κp i k B T ) (3a) P i LS→HS = 1 τ exp ( - D -k B T ln g k B T
) exp

( - E A + κp i k B T ) (3b)
where τ is a constant scaling factor, and κ a constant that relates the local pressure with the activation energy. We consider here that compressed springs produce a positive pressure (favouring the LS state), while elongated springs produces a negative pressure (favouring the HS state). After a Monte-Carlo standard procedure, one decides if a molecule switches or not. A Monte-Carlo electronic step (MCES) is completed when every molecule has been checked once, by comparing its transition probability (3a) or (3b) with a random number η ∈ (0, 1). A given molecule changes its state only if this number is smaller than the corresponding probability. After each MCES, the new equilibrium positions of all molecules are found.

In order to tune the different relative time scales of spin dynamics and lattice relaxation, we define the ratio

r = n 1 n 2
, where n 1 is the number of steps to solve differential equations ( 1), with probabilities (3) fixed and n 2 is the number of MCES per n 1 . It should be noted that a large r (large n 1 ) favours equilibrium distribution, allowing to accommodate to the distribution of local pressure p i , before any such pressure is changed, as previously considered in mechanoelastic model for reproducing static phenomena. It is in this way that the thermal hysteresis shown in figure 2, or slow relaxation curves 28 are simulated. Oppositely, a smaller r (high n 2 ) favours non-equilibrium distribution when the probabilities (3) change before the system accommodates to the distribution of local pressure, like in the case of fast phenomena subsequent to femtoseconds photoexcitation experiments 29 . This is due to the fact that molecular switching timescale (related to n 1 ) characterized by the timescale related to intra-molecular bond expansion (200 fs) is shorter than the typical lattice equilibration timescale (ns). In this paper, if not stated otherwise, we have used r = 1, that means that we solve one single step of differential equations after every MCES. This value corresponds to a situation when the lattice relaxation timescale is much longer than the molecular switching timescale.

For the present simulations, we use a rectangular shape sample, composed of 13982 molecules in triangular configurations which translates into real space in a 60 nm × 60 nm sample for the typical value of 1 nm per molecule. The value of the parameters considered in the simulations are in the range of those found by the estimation of the variation of the ligand field strengths in the two spin states due to the cooperativity and the bulk modulus 30 : κ = 0.2 × 10 -10 A m, k = 5 N m -1 , with the bond-length difference between the HS and LS molecules ∆r HL = 0.2 × 10 -10 m, while other parameters are µ = 0.7 N s m -1 , E A = 552 × 10 -23 J, ∆H = 1520 × 10 -23 J and g = 1096, thus giving a thermal transition temperature T 1/2 = ∆H/k B ln g of around 157 K (figure 2) close to the one experimentally reported for

[ F e III (3 -M eO -SalEen) ] P F 6 . We consider a temperature below the thermal transition temperature, at which most of the molecules are in the LS state. First, we run enough Monte Carlo steps in order to find the steady state at the corresponding temperature. Then the homogenous photoexcitation is simulated by choosing randomly a number of LS molecules, corresponding to the desired photoexcitation percentage, and by switching them to the HS state. Simultaneously, their individual temperature is increased with ∆T = 100 K 31 . This value is in line with the lattice heating experimentally observed of 15-20 K for a 15-20% photoexcitation. (consequently if 20% photoexcitation results in a global system heating of around 20 K, after the energy redistribution, we may assume that a 100% photoexcitation will result in a 100 K global heating, and, for practical purposes in Monte Carlo algorithm, we assume this value as individual increase of temperature for every photoexcited molecule). The temperature increase on the excited molecule allows to reduce the relaxation process towards its initial LS state due to lattice pressure. Here, as we are interested to simulate the behavior of the system for short time after photoexcitation when the elastic step appears and consider that the thermal diffusion is negligible and all molecules keep constant their temperatures attributed after the photoexcitation. This approximation was made for avoiding thermally-induced switching due to laser heating. We have also performed calculations taking into account the energy redistribution between the photoexcited molecules and their environment and, qualitatively, the same effects are observed. However, a correct description of the heat diffusion is more complex and this will be a further step towards more realistic simulations. Even if the heat diffusion could occur really quick for a high excitation (distance between two excited LS), the energy barriers between LS and HS are such that Arrhenius-like conversion takes much longer that the arrival of the elastic step 32 . Then, the system is allowed to evolve without any external influence by alternative checking molecules for switching and solving the system of differential equations.

The dependence of ∆X HS on the initial photoexcitation percentage is shown in figure 3 for a constant temperature T = 125 K. If the photoexcitation percentages are high enough, we notice a fast relaxation after the photoexcitation, followed by an increase of the HS fraction, which eventually surpasses the initial photoexcitation as observed experimentally. The variation of the HS fraction between the elastic step and after photoexcitation is ∆X HS (figure 1). However, if the initial photoexcitation percentage is too small, only the relaxation from HS to LS states can be observed and ∆X HS < 0. The apparent trough at ∆X HS = 0.15 is the manifestation of two competing and opposite going processes, namely the local relaxation and the non-local amplification.

The dependence of ∆X HS with temperature for a photoexcitation percentage of 15%, is shown in figure 4. If the temperature is close to the thermal transition temperature, then a large increase in the HS fraction is observed after photoexcitation. This increase corresponds to the experimentally observed elastic step 5 . The amplitude of this elastic step decreases with decreasing temperature until it vanishes at around 125 K. For temperatures below 125 K, only the relaxation subsequent to photoexcitation can be observed. By depicting the threshold at ∆X HS = 0, figure 4 (inset), we can conclude that with increasing the temperature, which stabilizes the higher entropy HS state, and thus the equilibrium HS fraction, the elastic step takes on a greater magnitude. This could be rationalized with the lattice getting softer, which requires less elastic cost for macroscopic volume expansion and molecular switching.

We plot in figure 5 the variation of the volume (that actually in our 2D case is just the area S) of the whole sample at 145 K for two opposite cases, showing or not the elastic step. As it takes some time for stress inside the lattice to propagate via an elastic distortion towards the edges of the system, the volume (area) variation does not exactly follow the X HS variation. In both cases, the initial HS fraction relaxation cannot be detected by the area variation: while some of the HS molecules relax, the area, which still corresponds to a state configuration before the photoexcitation, starts to increase in order to match the updated state due to internal pressure. Even in the case where only the relaxation is observed, the area still increases, as it should fit a state with ∆X HS higher than prior to photoexcitation. We should strengthen here that at short times described by our simulations only the propagation of phonons, driven by an harmonic-type potential influence the switch of molecules all around the sample, as the heating transfer between molecules is neglected and all of them keep their temperature. Only in case of high value of r, the area and the fraction would be somewhat artificially brought to concur in time. We emphasize that in our calculations r = 1, which corre- sponds to a realistic value that allowed reproducing all features on experimentally determined timescales.

The explanation of the experimental observation is strictly related to the microscopic evolutions inside the material, as it results from figure 6, where we present snapshots of the system taken at 145 K at different instants, before and during the elastic step. In figure 6 upper left, we show the situation just after photoexcitation, once molecules have switched from LS to HS state of higher volume but before the lattice relaxation starts. With this random distribution of HS molecules in a frozen lattice, most of them experience a very high pressure applied on them (blue circles). This high pressure will determine the switching back of several HS molecules to the lower volume LS state. Only the fact that the temperature of the photoexcited molecules has been increased prevents the switching back of all HS molecules (3). After a while (just at the end of relaxation part of X HS ), the relaxation of the lattice leads to a decreasing of pressure for molecules situated closer to the edge (red and black circles). This state is associated with the end of the relaxation part of the curves. Remarkably, the molecules still associated with a high pressure are situated inside a circular area, just in the middle of the sample. This observation shows actually how the relaxation starts from the edges, developing towards the center. Later on, during the elastic step, the lattice has almost completely relaxed. There are no more molecules with a high pressure acting on them, while the number of the molecules with a low pressure and in the HS state increased, especially towards edges and corners (black circles). This situation corresponds to figure 6.f, where one shows a "reorganized" sample: accumulations of HS molecules (small clusters) are concentrated towards corners and edges. Actually, this figure helps us to understand the nature of the elastic step. Several photoexcited molecules situated in the vicinity of corners / edge will lead to local deformations of the lattice resulting in creating more space for nearby LS molecules. This means that most springs around the LS molecules will be elongated and, therefore, the LS -HS switching probability for these molecules will increase, and many of them will switch to HS state. Therefore, the elastic step appears due to a reorganization of the molecules with different spin state and volume in the vicinity of corners / edges [movie in the Supplemental Material 33 ]. The same conclusions can be drawn from the representations of the maps of local pressures, i.e. compressions or elongations of every spring in the system, which are similar to local deformation fields in a continuous medium, as shown in Figure 7. Just after photoexcitation, the values of the local pressures are homogenously distributed all around the sample, with somewhat larger values around photoswitched molecules (Fig. 7,a); then high local pressures persists for a while in a central region of the sample, while they decrease faster for edges (fig. 7,b,c). Towards the end of the elastic step, the local pressures are again homogenous all around the sample, but with lower values.(Fig. 7,e,f).a,b,c,d,e,f correspond to different MCES shown in Figure 8.

The results presented din figure 8 supports these statements. They show, separately, the evolution of X HS for regions in the middle of the sample and towards the edge and support the hypothesis that most of the effect is due to the increase of the number of HS molecules towards the edge. While towards the edges there is almost no decrease of X HS after irradiation (as the additional pressure is rapidly freed by the increase of the volume), in the inside part the relaxation is important and the average of X HS is much smaller that for the whole sample. In the inset of figure 8, we present the variation of local pressures values during the elastic step. As qualitatively shown in figure 7, they decrease faster at the edges of the sample and slower in the bulk, but stabilize in the whole sample around the same non-null value towards the end of the elastic step.

The next question is to study the influence of the elastic constant on the magnitude of the elastic step. To clarify this point, we have simulated the behaviour of the system keeping constant all parameters of the system (as mentioned before), except the elastic constant. The temperature (145 K) and the photoexcitation percentage (20%) is also similar for all simulations. The results are presented in figure 9. At first sight, it can be stated that the elastic step exists even for smaller values of the elastic constant, when no hysteresis appears and that there is no monotonic evolution of the magnitude of the elastic step with the amplitude of the elastic constant (a maximum value of the elastic step is obtained for an intermediate elastic spring value). However, due to the different steady state before photoexcitation (see figure 9 inset), the interpretation of these results is rather difficult, as we have previously shown, the value of the elastic step depends also on the initial photoexcitation state. On one side, these results show the existence of an elastic step independently of the value of k, while we may imagine a chemical tuning of a material in order to obtain a larger elastic step. Different k entails different sound velocity, and indeed we observe that the onset of the elastic step in terms of MCES, not only its magnitude, is k dependent. The elastic step has a larger magnitude for an intermediate value of the elastic constant, as in the case of a soft lattice (low value of (k)), the local pressure is not large enough to cause the switching of many molecules, while in the case of a hard lattice (high value of (k))) the elastic distortions encounter more difficulties to propagate. [START_REF] Slimani | [END_REF] Next, we discuss the amplitude of the elastic step as a function of the ratio r between the number of steps in the differential equations solved, and the number of MCES between two successive solving of differential equation. As we have stated above, r is a measure for how fast is the lattice relaxation in comparison with individual molecular switching. The results for different r are presented in figure 10. If r is large then no relaxation can be observed after the photoexcitation, as the lattice relaxes (almost) completely till the next molecular switching, so the local pressure will diminish enough for all recently HS switched molecules and they will not switch back to LS. This situation corresponds to conventional LIESST experiments performed at low temperature where the photoinduced molecular HS state lives orders of magnitudes longer than the acoustic timescale. When r is somewhat smaller, the lattice does not completely relax and therefore some molecules (especially those situated in the middle of the lattice, according to figure 6) will still have a higher pressure and will switch back to LS state. This number of molecules is larger, so that is the relaxation, as r decreases more. However, in all situations, the FIG. 10: The influence of the ratio probability check / steps during molecular motion on the elastic step: the temperature is considered 145 K and a 20% photoexcitation. The final state is the same irrespective to the ratio, but the initial relaxation disappears if the lattice is let to relax after every molecular change.

value of the elastic step is similar -this can be easily understood if we take into account that it is the equilibrium state for a given spin configuration. The surface shows a faster variation for larger r values, while, for smaller values, an incubation time is observed, as experimentally reported.

A final point that we would like to address is the role of the homogeneity of the photoexcitation on the value of the elastic step. It was experimentally observed that in function of various parameters, like the absorption cross section of the material, the photoexcitation wavelength or sample size, the photoexcitation can be more or less homogenous 8 . We have therefore considered different initial photoexcited state configurations, from a complete inhomogeneous one, where only the molecules near one side are photoexcited (figure 11, top) to a complete homogeneous state, where the photoexcitation is completely random, as considered in previous paragraphs. Intermediate states with the number of photoexcited molecules decreasing exponentially from one side to another have also been investigated (figure 11 bottom). The last situation would be closest to the real situation of a single crystal, of which absorption obeys Lambert-Beer law [START_REF] Judd | Color in Business, Science and Industry[END_REF][START_REF] Enachescu | [END_REF] .

In figure 12 we show that the value of the elastic step obtained in the case of intermediate photoexcitation is three times higher than that in the case of homogeneous or inhomogeneous states, for similar global values of initial fraction of photoswitched molecules. This can be explained using the results presented in figure 10. If the photoexcitation is completely inhomogeneous and located on the surface, the photoexcitated molecules will only manifest in their immediate vicinity and will not have enough force to influence distant molecules. On the other side, in the case of completely random photoexcitation, most of the photoexcited molecules will be in the bulk part of the sample, so, either they will relax or wither they will determine the switching of a small number of LS. Therefore, the most favourable case is an intermediate one, as in figure 11 bottom. The proportion of molecules near the photoexcitation side is just high enough to determine the switching of almost all their nearby molecules, and the region till the photoexcitation population surpasses a threshold is large enough to guarantee their influence in a large part of the sample. As expected, the surface changes faster in the case of an inhomogenous state, as there are a higher number of photoexcited molecules near the edge, so the elastic distortion will need to cover a smaller distance, while for homogenous photoexcitation, the surface varies slowly.

III. CONCLUSIONS

In this paper, by using the mechanoelastic model, we have qualitatively reproduced the unusual self-amplified increase of the number of molecules in the high-spin state after femtosecond photoexcitation, recently reported in various experiments. The non-linear photo-response on acoustic timescale, for which elastic distortion propagates over the macroscopic sample, results from the competition between the local relaxation and the non-local amplification. We have explained this effect as a reorganisation of molecular states, resulting from the elastic distortion and due to the lattice deformations. Simulations points enhanced effect towards the edge and corners of the open boundary lattice. If the elastic lattice deformation origin of this self-amplification is well described by the present model, further work should consider size effects and heat FIG. 12: The evolution of nHS (up) and on the total surface (down) starting from different initial states with the same number of photoexcitated molecules but with different homogeneities (see figure 11). The maximum of the elastic step is attained for an intermediate non-homogenous initial state.

propagation, which manifest at a longer timescale, of mi-croseconds. A study of the dependence of the elastic step with the size will most likely lead to the conclusion that as the size of the crystals decreases the expansion occurs faster and elastic deformation is larger. Equally, a further study will concern the rigidity of the lattice, which may preclude amplification at low excitation density and favor it at high excitation density. Photothermal conversion is not discussed here. It depends on the proximity of the thermal phase transition, but it is limited by the thermal equilibrium timescale and energy redistribution from the photoexcited molecules to the global system should be taken into account. At least for systems homogeneously photoexcited, the description used in this paper remains valid but exponential decay of molecular temperature and exponential increase of the lattice temperature should be included. This will be a next step in our studies. Clearly one of the big hopes now is that we will be able to predict the elastic conversion time, whether it lags or not behind the volume (distortion) propagation, and that on the ultrafast time-scales for nanocrystals. We give away a taste for it in figure 11. The model can be improved by using parameters provided by ab-initio calculations such as vibration modes and frequencies or distortions during the switches between the states. We believe that the present model based on the coupling between the strain wave and the order parameter field is a new perspective extendable onto many other photoactive materials for which important molecular volume-change couple to lattice expansion [37][38][39][40][41][42] .

2 FIG. 1 :

 21 FIG.1: Experimental evidence of the elastic step. The time course of XHS following femtosecond laser excitation depends on the initial fraction of photoexcited molecules. For 2% (scaled ×2) of photoexcited molecules a relaxation occurs within 1 µs. For 7% a self-amplification occurs up to 27% within 1 ns during the elastic step.

FIG. 2 :

 2 FIG. 2: (main figure) Typical dynamical LS-HS transformation at 10000 MC steps/K; (inset) (up) LS molecule in LS environment and HS molecule in HS environment (down) The triangular network of particles connected with springs. In the initial state all the molecules are in the LS state (open circles). When a molecule switches its state to HS (red full circle), all the other molecules (small circles) change their positions. The lattice becomes then distorted: the distortions are bigger in the vicinity of switched molecule (short range effect), but they are present even in remote parts of the system (long range effect).

FIG. 3 :

 3 FIG. 3: Simulation of the evolution of the HS fraction after different percentages of photoexcitation showing the elastic step. Inset: the variations of ∆XHS after the elastic step as a function of photoexcitation ratio at constant T = 125 K.

FIG. 4 :

 4 FIG. 4: Simulation of the evolution of the HS fraction after photoexciting 15% of LS molecules showing the elastic step at different temperatures. Inset: the variations of ∆XHS after the elastic step as a function of temperature.

FIG. 5 :

 5 FIG. 5: Variation of XHS and of total area during elastic step at 145 K after two different initial photoexcitation percentages.

FIG. 6 :

 6 FIG. 6: Initial (upper left), intermediate and final states (lower right) after the propagation of the elastic distortionelastic step. Yellow molecules: low spin state. Other colours: high spin molecules with a colour scale corresponding to their local pressure: black: low pressure (very stable), red: moderate pressure (stable), blue: high pressure. The snapshots correspond to situations just after photoexcitation (a), during and after the relaxation (b, c) and during fluctuations around the elastic step (d, e, f), as shown in figure 8.

FIG. 8 :

 8 FIG. 8: Difference in the evolution of XHS for central and edge regions and comparison with the average XHS for the whole sample. Inset: variation of the local pressures for central and edge regions of the sample.

FIG. 11 :

 11 FIG. 11: Initial (left) and final state after the propagation of the elastic distortion-elastic step (right) for different initial states with 20% HS molecules. The same color scale as in figure 6 is used.

Acknowledgements This work was supported by National Research Council (CNCS-UEFISCDI), Romania, grant number PN-II-RU-TE-2014-4-0987, by the Institut Universitaire de France, Rennes Métropole, Région Bretagne (CREATE 4146), ANR (ANR-13-BS04-0002, CNRS (LIA IM-LED) and by Grants-in-Aid for Scientific Research C (No. 25400391, and No. 26400324) from MEXT of Japan.


Lucas Colabella

Adrián P Cisilino

Guillaume Häiat

Piotr Kowalczyk

MIMETIZATION OF THE ELASTIC PROPERTIES OF CANCELLOUS BONE VIA A PARAMETERIZED CELLULAR MATERIAL

) to mimic the elastic response of cancellous bone. Artificial microstructures are compared with actual bone samples in terms of elasticity matrices and their symmetry classes. The capability of the parameterized microstructure to combine the dominant isotropic, hexagonal, tetragonal and orthorhombic symmetry classes in the proportions present in the cancellous bone is shown. Based on this finding, two optimization approaches are devised to find the geometrical parameters of the artificial microstructure that better mimics the elastic response of a target natural bone specimen: a Sequential Quadratic Programming algorithm that minimizes the norm of the difference between the elasticity matrices, and a Pattern Search algorithm that minimizes the difference between the symmetry class decompositions. The pattern search approach is found to produce the best results. The performance of the method is demonstrated via analyses for 146 bone samples.

INTRODUCTION

Bones are hierarchical bio-composite materials with a complex multiscale structural geometry [START_REF] Carretta | Towards patient-specific material modeling of trabecular bone post-yield behavior[END_REF]). Bone tissue is arranged either in a compact pattern (cortical bone) or a spongy pattern (cancellous bone). Cancellous bone can be found in vertebral bodies and at the epiphyses of long bones. In the vertebral body, it is the main load bearing structure, where as in the appendicular skeleton, it transfers mechanical loads from the articular surface to cortical bone. Furthermore, trabecular bone quality is an important determinant of the overall bone strength and affects fracture risk. To better understand the mechanics of cancellous bone is of interest for the diagnosis of bone diseases (osteoporosis), the evaluation of the risk of fracture, and for the design of artificial bone [START_REF] Cowin | Bone Mechanics Handbook[END_REF].

Cancellous bone can be assimilated to a composite material with hierarchical structure. In a bottom-up description, the structure starts in the nanoscale (mineralized collagen fibril) and moves up to the sub-microscale (single lamella), the microscale (single trabecula), and mesoscale (trabecular bone) levels. Trabeculae are organized into a threedimensional lattice oriented mainly along the lines of stress, which forms a stiff and ductile structure that provides the framework for the soft bone marrow filling the intertrabecular spaces. Trabeculae consist of a nanometric extracellular matrix that incorporates hydroxyapatite, the bone mineral that provides bones rigidity [START_REF] Sansalone | Determination of the heterogeneous anisotropic elastic properties of human femoral bone: From nanoscopic to organ scale[END_REF][START_REF] Sansalone | Anatomical distribution of the degree of mineralization of bone tissue in human femoral neck: Impact on biomechanical properties[END_REF], and collagen, an elastic protein which improves fracture resistance [START_REF] Keaveny | Biomechanics of trabecular bone[END_REF].

Bone tissue mechanical properties and trabecular architecture are the main factors determining the mechanical properties of cancellous bone, which show a high dependency on species, anatomic site, age and size of the sample [START_REF] Fritsch | Universal" microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: Micromechanics-based prediction of anisotropic elasticity[END_REF][START_REF] Parkinson | Characterisation of Trabecular Bone Structure[END_REF]. The small dimensions of the trabeculae (of the order from tens to a cent of microns) hinder their mechanical characterization at tissuelevel. In recent years, nano-indentation has provided the means for the direct measurement of the elastic properties of trabecular bone tissue (a complete review of the available techniques, many of them indirect, can be found in the recent paper by [START_REF] Oftadeh | Biomechanics and Mechanobiology of Trabecular Bone: A Review[END_REF]). By means of high resolution nano-indentation, [START_REF] Brennan | Biomechanical properties across trabeculae from the proximal femur of normal and ovariectomised sheep[END_REF] studied the tissue property variations within a single trabecula; they found that Young´s modulus and hardness increase towards the core of the trabeculae. Despite this findings, it is the common assumption that the mechanical inhomogeneity and anisotropy of bone tissue has a minor impact on the apparent properties of cancellous bone and, consequently, it can be approximated by an isotropic tissue modulus (Kabel et al. 1999a, b).

Different experiments have shown that linear elasticity can predict the behavior of cancellous bone [START_REF] Keaveny | Trabecular bone exhibits fully linear elastic behavior and yields at low strains[END_REF]. The trabecular architecture determines the elastic anisotropy of cancellous bone, which can be described by a fourth rank elastic tensor, ℂ, which linearly relates stress and strains as 𝜎 = ℂ • 𝜀. The elastic tensor is determined in its most general form by 81 components. Cancellous bone is generally assumed to behave as an orthotropic structure in the mesoscale, with three planes of symmetry, what requires of only nine independent components to fully describe the elastic behavior of the structure [START_REF] Yang | The anisotropic Hooke's law for cancellous bone and wood[END_REF].

Relationships between elastic properties and structural parameters of cancellous bone have been proposed in the literature. Apparent density, 𝜌 𝑎𝑝𝑝 , which is the product of bone volume fraction, 𝐵𝑉/𝑇𝑉, and of bone tissue density, 𝜌, is the primary component affecting the mechanical properties of cancellous bone [START_REF] Oftadeh | Biomechanics and Mechanobiology of Trabecular Bone: A Review[END_REF]. Several studies have proposed linear and exponential forms that express the Young's modulus as a function of 𝜌 𝑎𝑝𝑝 , showing a high correlation with the experiments. However, these correlations exhibit limitations: the error in the elastic moduli can be up to 53% at certain volume fractions [START_REF] Van Rietbergen | Tissue stresses and strain in trabeculae of a canine proximal femur can be quantified from computer reconstructions[END_REF]) and they cannot individually predict the elastic properties of trabecular bone in different anatomic sites and species [START_REF] Helgason | Mathematical relationships between bone density and mechanical properties: A literature review[END_REF]. Similarly, cancellous bone elastic properties can be correlated to the fabric tensor, which is a descriptor of the anisotropy of trabecular bone [START_REF] Gross | Morphology-elasticity relationships using decreasing fabric information of human trabecular bone from three major anatomical locations[END_REF][START_REF] Maquer | Bone Volume Fraction and Fabric Anisotropy Are Better Determinants of Trabecular Bone Stiffness Than Other Morphological Variables[END_REF]. Another approach is to use large-scale finite element (FE) analysis of microstructural models built from micro-computed tomographic (CT) scans of real bone specimens. Finite element analyses can solve some of the drawbacks of the experimental techniques, since FE models can be subjected without restrictions to the load conditions needed to evaluate the anisotropic behavior of the microstructure. Moreover, the combination of finite element analysis (FEA) and modern in vivo high-resolution peripheral quantitative CT scanners, allow for the analysis of bone in-vivo in the peripheral skeleton [START_REF] Van Rietbergen | A survey of micro-finite element analysis for clinical assessment of bone strength: the first decade[END_REF]. FEA has been applied to large sets of data to find the orthotropic components of cancellous bone (Kabel et al. 1999a(Kabel et al. , 1999b)), which show that there are strong correlations between 𝐵𝑉/𝑇𝑉 and elastic and shear moduli, whereas this correlation is weak for the Poisson's ratio.

Another approach is to use parametric models of trabecular bone, which consist in artificial microstructures formed by plates and rods. Artificial microstructures may be criticized for being somewhat unrealistic. However, their main advantage is that the mesoscopic properties characterizing such microstructures can be expressed as explicit continuous functions of some well-defined geometrical parameters. Moreover, it has been found that models based on local morphometry, composed of individual rods and plates, help improving the understanding of local structural changes in the determination of bone stiffness (Stauber andMüller, 2006a, 2006b). Explicit relations between geometrical parameters and mesoscopic properties are crucial for modeling the microstructure evolution at the large scale -it allows to formulate the problem as merely the evolution of a set of scalar variables, which is much more efficient in terms of computational cost than the analysis of the geometric evolution of certain components of micro-CT-based actual bone microstructures. Examples in this sense are the artificial trabeculae developed by [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF] and [START_REF] Dagan | Single-trabecula building block for large-scale finite element models of cancellous bone[END_REF], which have been successfully employed in the modeling of long-term changes in morphological and mechanical properties of trabecular bone in the proximal femur, see Be'ery-Lipperman and Gefen ( 2005) and [START_REF] Kowalczyk | Simulation of orthotropic microstructure remodelling of cancellous bone[END_REF], respectively. Artificial trabecular microstructures have been also used for the developing of continuum models to describe the trabecular-bone stress-strain response [START_REF] Goda | A 3D elastic micropolar model of vertebral trabecular bone from lattice homogenization of the bone microstructure[END_REF]Goda and Ganghoffer 2015a) and multiaxial yield and failure criteria (Goda and Ganghoffer 2015b) by means of homogenization analyses. In turn, homogenized properties of parametric models of trabecular bone can be integrated into multiscale optimization methods for the design of bone substitutes and natural microscaffolds [START_REF] Hollister | Porous scaffold design for tissue engineering[END_REF] for tissue engineering. Such hierarchical structures can be fabricated by means of emerging 3D printing technologies [START_REF] Bose | Bone tissue engineering using 3D printing[END_REF][START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: A review[END_REF].

In this work the parameterized cancellous microstructures introduced by [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF] are analyzed in terms of their ability to mimic the elastic response of natural cancellous bone. Artificial microstructures are compared with actual bone samples in terms of their symmetry classes and their elasticity matrices represented in terms of the geometry parameters. Two optimization approaches are proposed in this paper to design the parameterized microstructure that better mimics the elastic response of a target natural bone specimen: a Sequential Quadratic Programming algorithm that minimizes the norm of the difference between the elasticity matrices, and a Pattern Search algorithm that minimizes the difference between the symmetry class decompositions. Both approaches use the geometry parameters as design variables.

ELASTIC PROPERTIES OF CANCELLOUS BONE

Experimental data

Two sources of experimental data are used: the database for human cancellous-bone specimens by Kabel et al. (1999aKabel et al. ( , 1999b)), and five bovine femoral bone specimens that are examined as part of this work.

The database by Kabel et al. (1999aKabel et al. ( , 1999b) ) provides the entire set of anisotropic elastic constants of 141 human cancellous-bone specimens of vertebral body, calcaneus, proximal tibia and distal femur. Specimen bone volume-to-total volume ratios cover the range 5% ≤ 𝐵𝑉/𝑇𝑉 ≤ 35%. The elastic constants are the results of finite element (FE) homogenization analyses performed on computer reconstructions of the specimen microarchitectures. Linear elastic and isotropic material properties were specified for the bone tissue, with a Young's modulus of 𝐸 = 1 𝐺𝑃𝑎 and Poisson's ratio 𝜈 = 0.3, so the homogenized results can be scaled for any value of the tissue modulus. The specimen imaging and homogenization procedures are fully described in Kabel et al. (1999b).

The five bovine samples were obtained from femoral bones. They were X-ray scanned using micro computed tomography (µCT) with a resolution of 17.7 𝜇𝑚. CT images were processed with BoneJ [START_REF] Doube | BoneJ: Free and extensible bone image analysis in ImageJ[END_REF]) to obtain geometrical data over the Regions of Interest (ROI). The results for bone volume fractions, trabecular thicknesses and trabecular spacings are reported in Table 1. The elastic modulus of the bone tissue was measured via microindentation tests using the [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF] method. Specimens were microindented in a TI 900 Triboindenter (Hysitron, MN, USA) using a Vickers diamond indenter. The maximum indentation load was set to 1500 𝑚𝑁, which was held constant for 45 𝑠 in order to minimize creep effects. The loading and unloading rates were set to 200 𝑚𝑁/𝑠 and 100 𝑚𝑁/𝑠 respectively. Eight indentations were performed on each sample. The resulting Young's modulus was 𝐸 𝑏 = 7.93 ± 0.86 𝐺𝑃𝑎. A Poisson's ratio 𝜈 = 0.3 was assumed. The elastic homogenization analyses were performed using the Fast Fourier Transform method with models built after the µCT geometry data. The procedure is described in [START_REF] Colabella | Calculation of Cancellous Bone Elastic Properties with the Polarization-based FFT Iterative Scheme[END_REF]. The resulting elasticity matrices 𝐶 𝑏 for the five specimens are reported in Appendix 1.

Elastic symmetry analyses

We propose here to study the elastic symmetries of the samples by decomposing their elasticity tensors into sums of orthogonal tensors belonging to the different symmetry classes. We use for this purpose the method by [START_REF] Browaeys | Decomposition of the elastic tensor and geophysical applications[END_REF]. This method relies on the following vectorial description of the elasticity tensor,

where 𝐶 𝑖𝑗 are the components of the elastic tensor ℂ in the Voigt notation. The normalization factors in the above expression are included so that the Euclidean norm of an arbitrary elastic tensor ℂ and its associated elastic vector 𝑿 are identical.

The vector description of the elastic tensor possesses the property that any symmetry class constitutes a subspace of a class of lower symmetry and an orthogonal projection on this subspace removes the lower symmetry part, see Table 2. Thus, when expressed in the so-called symmetry Cartesian coordinate system (see [START_REF] Cowin | On the identification of material symmetry for anisotropic elastic materials[END_REF], 𝑿 can be decomposed by a cascade of projections into a sum of vectors belonging to the symmetry classes triclinic, monoclinic, orthorhombic, tetragonal, hexagonal and isotropic:

𝑿 = 𝑿 𝑡𝑟𝑖 + 𝑿 𝑚𝑜𝑛 + 𝑿 𝑜𝑟𝑡 + 𝑿 𝑡𝑒𝑡 + 𝑿 ℎ𝑒𝑥 + 𝑿 𝑖𝑠𝑜 , (2)
The different elastic symmetry parts can be presented as fractions of the Euclidian norm of the elasticity vector ‖𝑿‖ as follows: Computations for the determination of the symmetry Cartesian coordinate system, the transformations into vector forms, the symmetry decompositions and normalizations were performed using the Matlab Seismic Anisotropy Toolkit (MSAT) by [START_REF] Walker | MSAT-A new toolkit for the analysis of elastic and seismic anisotropy[END_REF]. Figure 1 presents the results for the symmetry class decompositions of the 141 human bone samples in Kabel et al. (1999aKabel et al. ( , 1999b)). The cumulative decompositions in (4) are presented as functions of the sample 𝐵𝑉/𝑇𝑉. Table 3 summarizes the extreme values for the symmetry classes. It is observed from Figure 1 that the isotropic class accounts for the most significant fraction of the elasticity matrices over the complete 𝐵𝑉/𝑇𝑉 range.

𝑐 𝑖𝑠𝑜 =
Although its wide dispersion, the mean value of the isotropic fraction increases linearly with 𝐵𝑉/𝑇𝑉, from 𝑐 𝑖𝑠𝑜 ̅̅̅̅̅ = 0.49 for 𝐵𝑉/𝑇𝑉 = 5% to 𝑐 𝑖𝑠𝑜 ̅̅̅̅̅ = 0.69 for 𝐵𝑉/𝑇𝑉 = 35%. Its standard deviation is 𝑆𝐷 𝑖𝑠𝑜 = 0.12. The second relevant fraction is for the hexagonal class. Conversely to the isotropic class, the mean value of the hexagonal class decreases linearly with 𝐵𝑉/𝑇𝑉. The isotropic and hexagonal classes behave such that they add a constant, 𝑐 𝑖𝑠𝑜 ̅̅̅̅̅ + 𝑐 ℎ𝑒𝑥 ̅̅̅̅̅ =≅ 0.83, with a standard deviation 𝑆𝐷 𝑖𝑠𝑜+ℎ𝑒𝑥 = 0.07. The tetragonal class fraction is marginal; around 1% for nearly 98% of the samples. The orthorhombic class presents a wide dispersion, but its mean value is nearly constant 𝑐 𝑜𝑟𝑡 ̅̅̅̅̅ ≅ 0.10. The orthotropic symmetry,

𝑐 𝑜𝑟𝑡ℎ𝑜 = 𝑐 𝑜𝑟𝑡 + 𝑐 𝑡𝑒𝑡 + 𝑐 ℎ𝑒𝑥 + 𝑐 𝑖𝑠𝑜 , ( 5
)
presents a constant average value 𝑐 𝑜𝑟𝑡ℎ𝑜 ̅̅̅̅̅̅̅̅ = 0.93 with 𝑆𝐷 𝑜𝑟𝑡ℎ𝑜 = 0.04. This last result is consistent with the observation by [START_REF] Yang | The anisotropic Hooke's law for cancellous bone and wood[END_REF], who found that ℂ matrices present orthotropic symmetry with a 95% confidence level. by Kabel et al. (1999aKabel et al. ( , 1999b)) The results for the symmetry class decompositions of bovine samples are in Figure 2. Extreme values for the symmetry classes are in Table 3. It is observed that although extreme values of the bovine symmetry classes are within those of the human data ranges, the sum of the isotropic and hexagonal classes is 0.66 ≤ 𝑐 ℎ𝑒𝑥 + 𝑐 𝑖𝑠𝑜 ≤ 0.89, which has the lower limit outside the standard deviation of the human data. A similar behavior is observed for the orthotropic symmetry of the bovine samples, 0.79 ≤ 𝑐 𝑜𝑟𝑡ℎ𝑜 ≤ 0.91. The workflow in Figure 3 depicts the procedure for the development of the mimetic cancellous bone microstructure. The first step is the computation of the elasticity and BV/TV data of the parameterized cellular microarchitecture introduced by [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF], on which the development is based. The capabilities of the parameterized microstructure to mimic natural bone is assessed through the comparison and correlation of the elastic and BV/TV data with that of the natural specimens. Next, the elastic and BV/TV data are interpolated as functions of the geometrical parameters. These interpolations will play a key role in the implementation of the optimization algorithms in Section 4. Finally, the polynomial interpolations are checked for consistency. The details of this procedure are given next.

The parameterized cellular microstructure

The parameterized cellular microarchitecture introduced by [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF] is shown in Figure 4. It consists in a repeatable cell that is inscribed into a space-filling dodecahedron, so it can be arranged in rows and layers to completely fill the 3-D space. The geometry of the cell is described by Bezier curves and corresponding surface patches. Surface transitions between neighboring cells are smooth. Shaded areas denote trabecular surface while the hatched areas are the cross-sections at which the cell is "stuck" to identical neighboring cells.

The repeatable geometry is described in terms of four geometrical parameters: 𝑡 𝑐 , 𝑡 ℎ and 𝑡 𝑣 , which define proportions between trabecular plate widths and thicknesses to produce transversely isotropic microstructures in the 𝑥 1 -𝑥 2 plane; and 𝑡 𝑒 , which scales it in the 𝑥 1 direction to produce fully orthotropic microstructures. Parameters 𝑡 𝑐 , 𝑡 ℎ and 𝑡 𝑣 are non-dimensional as they are understood as fractions of the corresponding cell dimensions and may take values between 0 and 1. In order to produce feasible geometries, they must comply with the restrictions

𝑡 ℎ ≥ 𝑡 𝑐 , 𝑡 𝑣 ≥ 𝑡 𝑐 . ( 6
)
Parameter values can be set to produce microstructures with solid volume fractions in the range 0 < 𝐵𝑉 𝑇𝑉 ⁄ < 100%.

Homogenization analysis and comparison with natural specimens

Kowalczyk ( 2006) used FE homogenization to compute the elasticity matrices for a broad set of microstructures given by (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) quadruplets on the domains 𝑡 𝑐 ∈ (0,1), 𝑡 ℎ ∈ [𝑡 𝑐 , 1) and 𝑡 𝑣 ∈ [𝑡 𝑐 , 1) in increments of 0.05, 𝑡 𝑒 ∈ [0.6,1.4] in increments of 0.2, and elastic properties 𝐸 = 1 GPa and 𝜈 = 0.3. He studied in detail the dependency of the resultant effective elastic constants with 𝐵𝑉/𝑇𝑉, and tabulated functions for the elastic constants in terms of morphometric parameters: mean intercept length (MIL), volume orientation (VO) and star length distribution (SLD). From the comparison of the effective elastic constants to those of the natural samples by Kabel et al. (1999aKabel et al. ( , 1999b) ) he showed that individual ranges of 𝐸 𝑖𝑗 , 𝐺 𝑖𝑗 and 𝜈 𝑖𝑗 of the parameterized microstructures are always wider than those of the natural specimens for every 𝐵𝑉/𝑇𝑉 value.

In what follows, the above analysis is extended to assess the capability of the parameterized microstructure to mimic the natural elastic symmetries. To this end, the database of elastic constants for the parameterized microstructure was refined, such that the four geometric parameters were varied in increments of 0.05. The FE homogenization procedure of [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF] was used for the computation of the effective elastic constants. The repeatable cells that result for each combination of the geometric parameters were discretized with 8-node linear brick elements and appropriate boundary conditions that ensure fitting of all deformed neighboring cells to each other were specified. Six load cases were considered for each cell: pure stretching in three orthogonal directions 𝑥 1 , 𝑥 2 , 𝑥 3 (see Figure 4) and pure shear in three orthogonal planes (normal to the three directions), which were specified in terms of the displacement fields. Reaction forces were measured for each case and used to compute the elements of the elasticity matrix. Thus, the construction of the database consisted of 41,990 homogenization analysis that involved the solution of approximately 250,000 finiteelement models altogether.

Stiffness matrices of the parameterized microstructures are decomposed into their symmetry classes using the same procedure introduced earlier for the natural specimens.

The extreme values attained by the symmetry classes are reported in Table 3. It is found that the parameterized microstructure covers the complete extents of the tetrahedral and orthorhombic classes of the natural specimens. On the other hand, it fails to cover the lowermost values of the isotropic class, 0.24 ≤ 𝑐 𝑖𝑠𝑜 < 0.36, and the uppermost values of the hexagonal class, 0.49 > 𝑐 ℎ𝑒𝑥 ≥ 0.65 for the human specimens. However, it is worth noting that only a few samples lie within the excluded ranges: 5 samples have 𝑐 𝑖𝑠𝑜 < 0.36 and 2 samples have 𝑐 ℎ𝑒𝑥 > 0.49, i.e., less than 5% of the 141 samples in the database. Symmetry classes of the bovine samples lie always within the extents of the parameterized microstructure.

The capability of the parameterized microstructure to mimic the elastic behavior is further assessed in terms of the 𝐵𝑉/𝑇𝑉. Figure 5 shows the symmetry classes of the parameterized microstructures for the range of 𝐵𝑉/𝑇𝑉 of the natural samples; these are 17,522 data points (in order to keep the figure clear not all data points were plot). Results are presented for 𝑐 𝑖𝑠𝑜 , 𝑐 𝑖𝑠𝑜 + 𝑐 ℎ𝑒𝑥 and 𝑐 𝑜𝑟𝑡ℎ𝑜 in subfigures (a), (b) and (c), respectively. The corresponding data for the natural samples are also shown in the figures: gray areas indicate the standard deviation of the human samples (see Figure 1) while the square marks are the values of the bovine samples (see Figure 1). Figure 5(a) shows that with the only exception of the lowermost values, i.e., for 5% ≲ 𝐵𝑉 𝑇𝑉 ⁄ ≲ 7%, the parameterized microstructure is able to mimic the isotropic class of the natural trabeculae. Regarding 𝑐 𝑖𝑠𝑜 + 𝑐 ℎ𝑒𝑥 , Figure 5(b) shows that the parameterized microstructure completely encompasses the data of the natural samples (the gray area is hardly visible behind the symbols). Finally, the results for the orthotropic symmetry in Figure 5(c) show that, consistently with its geometric definition, the parameterized microstructure shows 𝑐 𝑜𝑟𝑡ℎ𝑜 = 1 over the complete 𝐵𝑉/𝑇𝑉 range, which results in a consistent overestimation of the orthotropic symmetry by the parameterized microstructure, the mean value and standard deviation of which is 6.6 ± 3.8% with respect to the human samples. Regarding the bovine samples, the overestimation ranges from 9% to 19%.

Polynomial interpolation

The discrete elastic-constant and symmetry-class data of the parameterized microstructures were examined to investigate their functionalities with the geometric parameters.

It was observed that coefficients of the stiffness tensor ℂ behave as continuous and smooth functions of 𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 and 𝑡 𝑒 over the complete range, and that in general, the 𝐶 𝑖𝑗 rise with the increment of the geometrical parameters. As examples, Figure 6 illustrates the behaviors of 𝐶 11 and 𝐶 12 as functions of 𝑡 𝑣 and 𝑡 𝑐 for 𝑡 ℎ = 0.6 and 𝑡 𝑒 = 1.2.

On the other hand, symmetry classes showed to be, in some cases, discontinuous functions of the geometric parameters. Figure 7 depicts the changes of the symmetry classes associated to the variation of the elasticity coefficients given in Figure 6. It can be observed that, although ℂ coefficients have a continuous variation, the hexagonal, tetragonal and orthorhombic symmetries present discontinuities.

Based on the above observations, the discrete elastic-constant data was used to interpolate an analytical expression for ℂ(𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ). The nine non-zero coefficients, 𝐶 11 , 𝐶 22 , 𝐶 33 , 𝐶 12 , 𝐶 13 , 𝐶 23 , 𝐶 44 , 𝐶 55 and 𝐶 66 , were interpolated polynomially by means of leastsquare fitting. Polynomials of order 5 to 12 were used; the quality of the interpolations was assessed in terms of the coefficient of determination (R2), the root-mean-square error (RMSE) and the residual sum of squares (SSres). Table 4 reports the averaged for the interpolations of the 9 coefficients. Since the 𝐶 𝑖𝑗 were interpolated separately, the thermodynamic requirement for ℂ(𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) was checked. The thermodynamic requirement (positive definiteness of strain energy) enforces the condition that the invariants of the elasticity matrix should be positive, or in other words, that both, ℂ(𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) and its inverse must be positive definite. These conditions were verified for the interpolated ℂ(𝑡 Bone volume-to-total volume ratio data was also interpolated polynomially to obtain an analytical expression for 𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ). Polynomials of order 2 to 5 were explored; the results are summarized in Table 5. It can be observed that the quality of the interpolation consistently improves with the polynomial order. The reliability of the interpolations was verified by checking that they result in positive 𝐵𝑉/𝑇𝑉 values. It was found that with the only exception of the interpolation of order 2, all the interpolations produce positive 𝐵𝑉/𝑇𝑉 values for every combination of the geometrical parameters within their validity intervals. The polynomial fitting of order 5 is selected for the rest of the work.

OPTIMIZATION

Figure 8: Workflow of the optimization analysis.

Problem statement

The optimization problem is to find the parameterized microstructure that better mimics the elastic response of a target natural bone specimen. The elastic equivalence among the microstructures is posed via two approaches: in terms of the elasticity matrix norm and in terms of its symmetry classes.

The workflow is illustrated in Figure 8. The analysis starts with the elastic homogenization of the target bone sample; the resultant elasticity matrix is decomposed into symmetry classes. The triclinic and monoclinic classes are suppressed to obtain the target elasticity matrix. Thus, the resultant ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 has orthotropic symmetry, and it is compatible with the artificial microstructure.

The problem posed in terms of the elasticity matrix norm consists in finding the quadruplet 𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 and 𝑡 𝑒 that minimizes the norm of the difference between the ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 and ℂ:

ℛ 1 = min ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 -ℂ(𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 )‖ ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖ , (7)
The problem posed in terms of the symmetry classes consists in finding the quadruplet 𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 and 𝑡 𝑒 that minimizes the overall difference between the symmetry class decompositions:

ℛ 2 = min √ (𝑐 𝑖𝑠𝑜 𝑡𝑎𝑟𝑔𝑒𝑡 -𝑐 𝑖𝑠𝑜 ) 2 + (𝑐 ℎ𝑒𝑥 𝑡𝑎𝑟𝑔𝑒𝑡 -𝑐 ℎ𝑒𝑥 ) 2 +(𝑐 𝑡𝑒𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 -𝑐 𝑡𝑒𝑡 ) 2 + (𝑐 𝑜𝑟𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 -𝑐 𝑜𝑟𝑡 ) 2 , (8)
where 𝑐 𝑖𝑠𝑜 𝑡𝑎𝑟𝑔𝑒𝑡 , 𝑐 ℎ𝑒𝑥 𝑡𝑎𝑟𝑔𝑒𝑡 , 𝑐 𝑡𝑒𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 , 𝑐 𝑜𝑟𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 and 𝑐 𝑖𝑠𝑜 , 𝑐 ℎ𝑒𝑥 , 𝑐 𝑡𝑒𝑡 , 𝑐 𝑜𝑟𝑡 are the normalized symmetry classes of the target and the parameterized microstructures, respectively; 𝑐 𝑖𝑠𝑜 , 𝑐 ℎ𝑒𝑥 , 𝑐 𝑡𝑒𝑡 and 𝑐 𝑜𝑟𝑡 are functions of 𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 and 𝑡 𝑒 .

For the two approaches, geometric parameters are subjected to the inequality restrictions in equation ( 6) and at the same time, they must comply with the restriction imposed by the bone volume-to-total volume ratio of the target microstructure,

𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) = 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 , (9)
Alternatively, the restriction of the bone volume-to-total volume ratio might be relaxed, such as

𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 (1 -𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 ) ≤ 𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) ≤ 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 (1 + 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 ), ( 10
)
where 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 is a prescribed tolerance that could take any value in the interval [0,1].

Algorithms

Based on the continuity of the interpolations of ℂ and its symmetry classes (see Section 3.3), different algorithms are used for the optimization problems posed in equations ( 7) and (8).

Being the interpolation of the ℂ coefficient continuous, the optimization problem ( 7) is solved using and active-set Sequential Quadratic Programming (SQP) algorithm. SQP methods solve a sequence of optimization subproblems, each of which optimizes a quadratic model of the objective function subjected to a linearization of the constraints. We use in this work the SQP algorithm of the Matlab Optimization Toolbox (fmincon), which is a customized implementation of the algorithm by [START_REF] Gill | Procedures for optimization problems with a mixture of bounds and general linear constraints[END_REF][START_REF] Gill | Numerical Linear Algebra and Optimization[END_REF]. The constrain functions are used for the restriction on 𝐵𝑉/𝑇𝑉. In the case of a prescribed target 𝐵𝑉/𝑇𝑉, the equality function constrain, 𝑐 𝑒𝑞 (𝑥) = 0, is

𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) -𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 = 0, ( 12
)
while if 𝐵𝑉/𝑇𝑉 is defined as in equation ( 10), two 𝑐 𝑖𝑛𝑒𝑞 (𝑥) ≤ 0 are specified:

𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) -𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 (1 + 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 ) ≤ 0 (13)
and

𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 (1 -𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 ) -𝐵𝑉 𝑇𝑉 ⁄ (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) ≤ 0. ( 14
)
Inequality constrains are the restrictions on the geometric parameters in equation ( 6). Thus, 𝐴 𝑖𝑛𝑒𝑞 (𝑥) • 𝑥 ≤ 𝑏 𝑖𝑛𝑒𝑞 is

[ 1 -1 0 0 1 0 -1 0 ] [ 𝑡 𝑐 𝑡 ℎ 𝑡 𝑣 𝑡 𝑒 ] ≤ [ 0 0 ]. ( 15
)
Finally, inequality constrains 𝑏 𝑙𝑜𝑤 ≤ 𝑥 ≤ 𝑏 𝑢𝑝 are used to set the validity intervals for the polynomial interpolations in Section 3.2. For the polynomials of order 10, these are Stopping criteria consist of two tolerances: FunctionTolerance, a lower bound on the change in the value of the objective function during a step, and StepTolerance, the termination tolerance on the step size. Iterations ends when either of the conditions is achieved. Both tolerances were set with the default values of 10 -6 .

𝑏 𝑙𝑜𝑤 = [ 0.
The optimization posed in term of the symmetry classes is solved using the derivativefree constrained direct search solver patternsearch (PS) of the Matlab Global Optimization Toolbox. patternsearch computes a sequence of points that approach an optimal. At each step, the algorithm searches a set of points, called a mesh, around the current point. The mesh is formed by adding the current point to a scalar multiple of a set of vectors called a pattern. If the pattern search algorithm finds a point in the mesh that improves the objective function at the current point, the new point becomes the current point at the next step of the algorithm.

The problem for the PS algorithm is specified using the same outline for the SQP in equation ( 11), but with the objective function 𝑓(𝑥) = ℛ 2 (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) of equation ( 8).

Restrictions to the geometric parameters due to the feasibility of the microstructures, 𝐵𝑉/𝑇𝑉 and confidence intervals for the polynomial interpolations are the same as for the SQP algorithm. Stopping criteria involves tolerances for FunctionTolerance, the difference between the function value at the previous best point and function value at the current best point; MeshTolerance, the minimum size for the search mesh; and

StepTolerance, the minimum distance from the previous best point to the current best point. The three tolerances were set with the default values of 10 -6 .

Verification and tuning

The optimization procedures were verified, tested and tuned by assessing its effectiveness to identify microstructures among those of the database used for the elasticity-matrix polynomial fitting (see Section 3.3). To this end, 100 parameterized microstructures were randomly selected from the database to serve as target microstructures. The optimization problems were solved with the SQP and PS methods for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% and 5%. Since the thermodynamic requirement was checked only for discrete combinations (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) in Section 3.3, there is no guarantee that all possible combinations will satisfy it. Therefore, the thermodynamic requirement for ℂ(𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) was checked at the end of each optimization procedure.

Preliminary tests had shown that the performances of both algorithms were sensitive to the initial values (seeds) of the geometric parameters. Thus, SQP optimizations were attempted with different seeds as many times as necessary until the objective functions attained the condition ℛ 1 < 0.001; the PS optimizations were run four times for different sets of random seeds and the best result reported.

𝑒 𝐵𝑉 𝑇𝑉 ⁄ = 𝐵𝑉 𝑇𝑉 ⁄ -𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑎𝑟𝑔𝑒𝑡 , ( 19
)
are 𝑒 𝐵𝑉 𝑇𝑉 ⁄ ̅̅̅̅̅̅̅̅̅ = 6 • 10 -5 and -5 • 10 -4 for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% and 5%, respectively.

The PS optimization produces accurate results for the symmetry classes, with maximum errors of a few percent for 𝑒 𝑡𝑒𝑡 . Error 𝑒 𝑡𝑒𝑡 diminishes from 8% to 3% with the increment of 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 . However, this is not the case for all the symmetry classes; note that 𝑒 ℎ𝑒𝑥 and 𝑒 𝑜𝑟𝑡 present their minima for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1%. Maximum standard deviations are also for 𝑒 𝑡𝑒𝑡 , and they reduce from 50% to 17% with the increment of 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 . The largest error in the geometric parameters is for 𝑡 𝑐 , which diminishes from 17% to 5% as the tolerance for the bone volume fraction is relaxed from 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 0 to 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 5%. In contrast, minimum errors for 𝑡 ℎ , 𝑡 𝑣 and 𝑡 𝑒 are for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 0, and they deteriorate with 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 ; in any case, maximum errors are of a few percent. The standard deviations for the geometric parameters are not sensitive to 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 and they are from 25% to 40%. Like for the SQP procedure, mean 𝐵𝑉 𝑇𝑉 ⁄ of the optimal microstructures are very close the target values, in these cases 𝑒 𝐵𝑉 𝑇𝑉 ⁄ ̅̅̅̅̅̅̅̅̅ = 7 • 10 -4 and -5 • 10 -3 for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% and 5%, respectively.

In what respects to the computational performance, SQP algorithm had to be run, in average, 1.4 times for different seeds to attain the goal ℛ 1 < 0.001. Each run comprised an average of 32 iterations and 180 function evaluations. The performance of the SQP algorithm was independent of 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 . The PS algorithm employed on average 3.7 iterations with 1720 function evaluations for each of the seed value sets. It is worth noting that, although the larger number of function evaluations, the PS method was always faster than the SQP. This is because, besides the greater complexity of the algorithm, the SPQ method needs of the cost function gradient, which is around five times more expensive to evaluate than the evaluation of the cost function itself.

The above results show that SQP and PS methods are effective to retrieve randomly selected specimens from the parameterized microstructure database. The SQP method produces accurate results in terms of the geometric parameters, but they are poor in what respects to the symmetry classes. The PS method produces accurate results in terms of the symmetry classes, which, at the same time, show good results for the geometric parameters. The two methods are further assessed next, by finding parameterized microstructures that better mimic the elastic response of natural microstructures.

APPLICATION TO NATURAL TRABECULAR SAMPLES

Human specimens

Elasticity matrices of the 141 human bone samples were filtered to retrieve their orthotropic parts. To do so, components 𝐶 14 , 𝐶 15 , 𝐶 16 , 𝐶 24 , 𝐶 25 , 𝐶 26 , 𝐶 34 , 𝐶 35 , 𝐶 36 , 𝐶 45 , 𝐶 46 and 𝐶 56 were set equal to zero for the elastic tensors oriented in the symmetry Cartesian coordinate system and their symmetry class decompositions were computed such that 𝑐 𝑜𝑟𝑡 𝑜𝑏𝑗 + 𝑐 𝑡𝑒𝑡 𝑜𝑏𝑗 + 𝑐 ℎ𝑒𝑥 𝑜𝑏𝑗 + 𝑐 𝑖𝑠𝑜 𝑜𝑏𝑗 = 1. The resultant data was used as target values for the SQP and PS optimization procedures. Optimizations were performed for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 0%, 1%, 5% and 10%.

Symmetry classes for the obtained effective elastic properties are assessed in relation to the mean values of their extreme fractions for the human microstructures in Table 3, this is:

𝑒 𝑐 𝑖𝑠𝑜 , = 𝑐 𝑖𝑠𝑜 -𝑐 𝑖𝑠𝑜 𝑡𝑎𝑟𝑔𝑒𝑡 1 2 (𝑐 𝑖𝑠𝑜 𝑚𝑎𝑥 + 𝑐 𝑖𝑠𝑜 𝑚𝑖𝑛 ) , 𝑒 𝑐 ℎ𝑒𝑥 , = 𝑐 ℎ𝑒𝑥 -𝑐 ℎ𝑒𝑥 𝑡𝑎𝑟𝑔𝑒𝑡 1 2 (𝑐 ℎ𝑒𝑥 𝑚𝑎𝑥 + 𝑐 ℎ𝑒𝑥 𝑚𝑖𝑛 ) , 𝑒 𝑐 𝑡𝑒𝑡 , = 𝑐 𝑡𝑒𝑡 -𝑐 𝑡𝑒𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 1 2 (𝑐 𝑡𝑒𝑡 𝑚𝑎𝑥 + 𝑐 𝑡𝑒𝑡 𝑚𝑖𝑛 ) and 𝑒 𝑐 𝑜𝑟𝑡 , = 𝑐 𝑜𝑟𝑡 -𝑐 𝑜𝑟𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 1 2 (𝑐 𝑜𝑟𝑡 𝑚𝑎𝑥 + 𝑐 𝑜𝑟𝑡 𝑚𝑖𝑛 ) . ( 20
)
Neither the SQP nor the PS showed significant improvements with the relaxation of the 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 constrain, so only the results for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% are presented next. Mean errors for 141 samples when using the SQP method are 𝑒 𝑐 𝑖𝑠𝑜 , ̅̅̅̅̅ = 0.23, 𝑒 𝑐 ℎ𝑒𝑥 , ̅̅̅̅̅̅ = -0.50, 𝑒 𝑐 𝑡𝑒𝑡 , ̅̅̅̅̅̅ = 0.54 and 𝑒 𝑐 𝑜𝑟𝑡 , ̅̅̅̅̅̅ = 0.17, while for the PS method they are 𝑒 𝑐 𝑖𝑠𝑜 , ̅̅̅̅̅ = 0.07, 𝑒 𝑐 ℎ𝑒𝑥 , ̅̅̅̅̅̅ = -0.09, 𝑒 𝑐 𝑡𝑒𝑡 , ̅̅̅̅̅̅ = -0.03 and 𝑒 𝑐 𝑜𝑟𝑡 , ̅̅̅̅̅̅ = -0.02. It can be observed that, like in Section 4.3, the PS achieved better results than the SQP. Figure 9 presents the results for the PS optimizations in terms of the relative error for the elasticity matrices,

𝑒 ‖ℂ‖ = ‖ℂ‖ -‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖ ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖ . ( 21
)
It can be observed that 𝑒 ‖ℂ‖ ̅̅̅̅̅ and its dispersion diminish with 𝐵𝑉/𝑇𝑉, from 0.3 ≲ 𝑒 ‖ℂ‖ ≲ 4 for 𝐵𝑉/𝑇𝑉 < 10% to 0.2 ≲ 𝑒 ‖ℂ‖ ≲ 0.8 for 𝐵𝑉/𝑇𝑉 > 32%. At the same time, it is interesting to note that, with only a few exceptions, 𝑒 ‖ℂ‖ > 0, what implies that the optimized parameterized microstructures are, in general, stiffer than the target natural ones. It might be argued that the parameterized microstructures make, in terms of stiffness, a more efficient use of the material than natural microstructures, being this greater efficiency more noticeable for low 𝐵𝑉 𝑇𝑉 ⁄ . This behavior allows to explain the poor performance of the SQP optimization: since parameterized microstructures are in general stiffer than the natural microstructures for a given 𝐵𝑉 𝑇𝑉 ⁄ , the objective function posed in terms of the elastic matrices (see ( 7)) is hard to minimize. On the other hand, the PS optimization, which is posed in terms of the symmetry classes, is not affected by the overall stiffness, and thus it conducts to better results. The 𝐵𝑉/𝑇𝑉 results confirm this observation: mean values of the error in ( 19) are close to zero for the PS method irrespectively of 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 (𝑒 𝐵𝑉 𝑇𝑉 ⁄ ̅̅̅̅̅̅̅̅̅ = -1 • 10 -3 , 5 • 10 -3 and 2 • 10 -2 for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% , 5% and 10% , respectively), whereas for the SQP they are always almost coincident with the lowest value of the tolerance range. Further results concerning the performance of the SQP and PS algorithms are given in the Appendix.

Bovine specimens

The PS optimization was performed for the orthotropic part of the stiffness matrices of the 5 bovine specimens described in Section 2.1. Tolerance for the volume fraction was set 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1%.

The optimized microstructures are shown in Figure 10 together with their corresponding target natural samples. The stiffness matrices of the parameterized microstructures are given in the Appendix; their geometric parameters and 𝐵𝑉 𝑇𝑉 ⁄ values are reported in Table 8. Table 9 and Table 10 report the errors for the symmetry classes and the elasticity coefficients, respectively. Error for the symmetry classes are in relation to the ranges of the extreme fractions for the bovine microstructures in Table 3 as in equation ( 20).

Table 9 shows that, with the only exception of the tetragonal symmetry, errors for the symmetry classes are very low. The large relative errors for the tetragonal symmetry are explained due to its small relative contribution to the elasticity matrix (see Figure 2). The worst performance is for Sample #2, which has the particularity of having the lowest isotropic class and the highest orthorhombic class fractions, 𝑐 𝑖𝑠𝑜 = 0.54 and 𝑐 𝑜𝑟𝑡 = 0.27 , respectively. Bone volume fractions of the mimetic microstructures are coincident to those of the natural specimens (see Table 1) up to the second significant figure, which is coherent to the 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 1% used for the analyses.

Error for the elasticity coefficients in Table 10 were computed after the scaling of the Young´s modulus to make 𝑒 ‖ℂ‖ = 0. The resultant scaling factors range 0.41 < 𝐸′ 𝐸 < 0.72 ⁄ with a mean value of 𝐸′ 𝐸 ⁄ ̅̅̅̅̅̅̅ = 0.55, which coincides with that of the human sample analysis in the previous section. It is observed that as for the human samples, 𝐶 12 presents the highest mean error level (around 9%); maximum errors are of around 17% for 𝐶 11 and 𝐶 22 of Sample #5.

Sample

Geometric

CONCLUSIONS

This work introduces a procedure for the design of artificial parameterized microstructures that mimic the elastic response of cancellous bone. The procedure is based on the parameterized microstructure by [START_REF] Kowalczyk | Orthotropic properties of cancellous bone modelled as parameterized cellular material[END_REF], the geometric parameters of which are optimized to minimize the differences between the symmetry classes of the target and the artificial microstructure elastic tensor.

Symmetry class analyses of experimental data from Kabel et al. (1999aKabel et al. ( , 1999b) ) and of specimens processed as part of this work show that elastic symmetries can be related to the specimen 𝐵𝑉/𝑇𝑉. The isotropic symmetry class constitutes the main fraction of the elastic tensor; it increases linearly with the specimen bone volume fraction, from around 50% for 𝐵𝑉/𝑇𝑉 = 5% to 70% for 𝐵𝑉/𝑇𝑉 = 35%. The isotropic and hexagonal classes add to a constant, such that they account for around 82% of the elastic tensor over the complete 𝐵𝑉/𝑇𝑉 range. The orthotropic symmetry, given by the addition of the isotropic, hexagonal, tetragonal and orthorhombic classes, constitute around 93% of the elastic tensor, independently of the 𝐵𝑉/𝑇𝑉.

The parameterized artificial microstructure is orthotropic by construction. It is shown in this work that it has the capability to combine the isotropic, hexagonal, tetragonal and orthorhombic symmetry classes in the proportions present in the cancellous bone. Analytical expressions for the elastic matrix in terms of the microstructure geometrical parameters are provided. These expressions could be integrated into multiscale design methodologies and used to explore the design of bone substitutes and natural microscaffolds. Free Material Optimization Methods (FMO) are seen as a promising approach in this sense.

Two optimization methods are proposed to find the parameterized microstructure that better mimics the elastic response of a target natural bone specimen: a Sequential Quadratic Programming algorithm to minimize the difference between the elasticity matrices, and a Pattern Search algorithm to minimize the difference between the symmetry class decompositions. Both approaches use the geometry parameters as design variables, polynomial interpolations to evaluate the parameterized microstructure elasticity matrix, and 𝐵𝑉/𝑇𝑉as a restriction. The Pattern Search approach is found to produce the best results. The analyses of 146 natural cancellous bone specimens resulted in mimetic microstructures whose symmetry class decompositions differ on average 6% with respect to the target values.

The results for the elasticity matrix error allows to observe that the optimized microstructures are in general stiffer than their natural counterparts; this behavior is more noticeable for low 𝐵𝑉/𝑇𝑉. This deviation can be compensated by selecting the Young's modulus for the optimized microstructure material such that norm of the difference between elasticity matrices of the target and optimized microstructure vanishes. Clearly, the Young´s modulus scaling does not affect the elastic symmetries. The mean value for such scaling factor was found equal to 0.55, i.e., the parameterized microstructure material should have, in average, half the stiffness of the trabecular bone tissue. After scaling, average errors between the optimized and target elasticity matrix coefficients do not exceed 10% relative to the matrix norm.

Additional research is needed to further assess the effectivity of the parameterized microstructures to mimic the behavior of cancellous bone elastic response. A promising approach is to compare the elastodynamic responses of natural and parameterized microstructures by means of ultrasound analyses, and to correlate their behaviors with the elastic symmetry classes and the geometrical parameters. In turn, these results could be used to refine the optimization criteria of the Pattern Search method. The PS optimization resulted in a mean value for the residual ℛ 2 ̅̅̅̅ = 0.06 with standard deviation 𝑆𝐷 ℛ 2 = 0.09. Like for the SQP approach, the relaxation of the 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 constrain did not result in significant improvements for ℛ 2 . In addition to the error for the elasticity matrices shown in Figure 9, Figure 12 presents the errors (18) for the symmetry classes. Like for the SQP optimization, 𝑒 𝑐 𝑡𝑒𝑡 presents the largest dispersion and there are tendencies to overestimate 𝑐 𝑖𝑠𝑜 and to underestimate 𝑐 ℎ𝑒𝑥 . The high dispersion of 𝑒 𝑐 𝑡𝑒𝑡 is consequence of the small relative contribution of the tetragonal symmetry to the elasticity matrix, which is always 𝑐 𝑡𝑒𝑡 < 0.05, see Table 3 and Figure 1.

Bovine Samples

Target elastic matrices of the bovine femoral samples The PS optimization of each sample was run ten times using different seeds. For each case, the best two outcomes (these are, the solutions that achieve the lowest values of the objective function) were compared to assess the repeatability of the solutions. It was found that for Samples #2 and #4, the residuals of the best two solutions are almost coincident (they differ less than 1%), while the resultant values for the geometrical parameters coincide within 2%. For samples #1 and #3 the lowest two residuals differ in around 12% and the geometrical parameters have discrepancies of up to 80%. It is interesting to note that the parameter 𝑡 𝑒 , the one which governs microstructure orthotropy, showed a remarkable repeatability, it presented discrepancies within 0.2% for the analyses of Samples #1, #2, #3 and #4. In contrast, for Sample #5, residuals of the best two solutions differ in nearly 90% and the geometrical parameters up to 70%; the best performance is for 𝑡 𝑒 , which presents a discrepancy discrepancy of 13%. The behavior for 𝑡 𝑒 could be explored as a mean for the refinement of the optimization procedure.

Finally, the resultant elastic matrices for the mimetic parameterized microstructures are:

Figure 1 :

 1 Figure 1: Symmetry class decomposition of the elasticity matrices of the 141 human-bone specimens reportedbyKabel et al. (1999aKabel et al. ( , 1999b)). Error bars indicate the standard deviations from the interpolated mean values.

Figure 2 :Figure 3 :

 23 Figure 2: Symmetry class decomposition of the bovine samples.

Figure 4 :

 4 Figure 4: Geometry of a repeatable cell.

Figure 5 :

 5 Figure 5: Comparison of elastic symmetry classes of the natural and the parameterized trabecular microstructures: (a) isotropic class, (b) isotropic + hexagonal classes, (c) orthotropy.

Figure 6 :

 6 Figure 6: Parameterized-microstructure elasticity coefficients 𝐶 11 and 𝐶 12 as functions of 𝑡 𝑐 and 𝑡 𝑣 for 𝑡 ℎ = 0.6 and 𝑡 𝑒 = 1.2.

Figure 7 :

 7 Figure 7: Parameterized-microstructure symmetry classes as functions of 𝑡 𝑣 and 𝑡 𝑐 for 𝑡 ℎ = 0.6 and 𝑡 𝑒 = 1.2.

Figure 9 :

 9 Figure 9: PS optimization of the human samples: stiffness-matrix error as function of 𝐵𝑉 ⁄ 𝑇𝑉. The solution of the PS optimization can be improved for 𝑒 ‖ℂ‖ by doing a convenient selection of the Young's modulus of the artificial microstructure material. If the artificial microstructure is constructed using a material with Young´s modulus 𝐸 ′ = ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖ ‖ℂ‖ 𝐸 , the error for the elastic matrix norm 𝑒 ‖ℂ‖ = 0. Clearly, this scaling of the Young´s modulus does not affect the elastic symmetries. The above analysis was performed for the 141 samples. The average scaling factor for the Young´s modulus was found 𝐸′ 𝐸 ⁄ ̅̅̅̅̅̅̅ = 0.55 with and standard deviation 𝑆𝐷 𝐸′ 𝐸 ⁄ = 0.20. Resultant mean errors and standard deviations for

Figure 10 :

 10 Figure 10: Natural bovine specimens and their mimetic parameterized microstructures.

Figure 11 :

 11 Figure 11: SQP optimization of human samples: (a) residuals and (b) symmetry-class errors as functions of 𝐵𝑉 ⁄ 𝑇𝑉.

  Figure 12 allows to observe that mean errors for the symmetry classes are much lower than those of the SQP approach in Figure 11(b) : 𝑒 𝑐 𝑖𝑠𝑜 ̅̅̅̅̅ = 0.06, 𝑒 𝑐 ℎ𝑒𝑥 ̅̅̅̅̅̅ = -0.13, 𝑒 𝑐 𝑡𝑒𝑡 ̅̅̅̅̅̅ = -0.06 and 𝑒 𝑐 𝑜𝑟𝑡 ̅̅̅̅̅̅ = -0.02; standard deviations are 𝑆𝐷 𝑒 𝑖𝑠𝑜 = 0.11, 𝑆𝐷 𝑒 ℎ𝑒𝑥 = 0.27, 𝑆𝐷 𝑒 ℎ𝑒𝑥 = 0.27 and 𝑆𝐷 𝑜𝑟𝑡 = 0.18.

Figure 12 :

 12 Figure 12: PS optimization of human samples: symmetry-class error as function of 𝐵𝑉 ⁄ 𝑇𝑉.

Table 1 : Geometrical data of the bovine samples

 1

	Specimen	ROI	Pixel	BV/TV	Trabecular	Trabecular	Normalized
		dimensions	size	[%]	thickness,	spacing,	trabecular
		[𝑚𝑚]	[𝜇𝑚]		𝑡 [𝜇𝑚]	𝑠 [𝜇𝑚]	thickness, 𝑡/(𝑡 +
							𝑠)
	1	8.7x9.6x9.5	17.7	25	157	491	0.24
	2	9.9x10.3x9.8	17.7	38	241	471	0.34
	3	9.8x10.1x7.7	17.7	30	180	621	0.22
	4	9.6x9.2x8.9	17.7	20	180	757	0.19
	5	9.7x9.5x7.6	17.7	21	145	573	0.20

Table 2 :

 2 Number of planes of symmetry and dimension of the subspaces. Note that the dimension number corresponds to the number of distinct components of the elastic tensor.

	Symmetry class Planes of symmetry Dimension
	Triclinic	0	21
	Monoclinic	1	13
	Orthorhombic	3	9
	Tetragonal	5	6
	Hexagonal	7	5
	Isotropic	∞	2

Table 3 :

 3 Extreme values of the symmetry classes of the natural and parameterized trabecular microstructures.

				Bovine Samples	Parameterized
				(this work)	Kowalczyk (2006)
		Min	Max	Min	Max	Min	Max
	𝑐 𝑖𝑠𝑜	0.24	0.84	0.50	0.79	0.36	1.00
	𝑐 ℎ𝑒𝑥	0.02	0.65	0.04	0.20	0.00	0.49
	𝑐 𝑖𝑠𝑜 + 𝑐 ℎ𝑒𝑥	0.56	0.95	0.66	0.89	0.46	1.00
	𝑐 𝑡𝑒𝑡	0.00	0.05	0.00	0.01	0.00	0.15
	𝑐 𝑜𝑟𝑡	0.00	0.33	0.02	0.13	0.00	0.50
	𝑐 𝑜𝑟𝑡ℎ𝑜	0.81	0.99	0.79	0.91	1.00	1.00
	𝑐 𝑚𝑜𝑛	0.00	0.08	0.00	0.04	-	-
	𝑐 𝑡𝑟𝑖	0.01	0.19	0.08	0.19	-	-

. Error bars indicate the standard deviations from the interpolated mean values.

Symmetry class Human Samples

Kabel et al (1999)

Table 4 :

 4 Quality assessment of the polynomial interpolations of the elasticity coefficients in terms of the coefficient of determination (R

2 ), the root-mean-square error (RMSE) and the residual sum of squares (SSres).

  𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ) for quadruplets in the intervals 𝑡 𝑐 ∈ [0.05, 0.95], 𝑡 ℎ ∈ [𝑡 𝑐 , 0.95], 𝑡 𝑣 ∈ [𝑡 𝑐 , 0.95] and 𝑡 𝑒 ∈ [0.6,1.4] in increments of 0.01. The thermodynamic requirement was found valid on most of the interpolation range; validity intervals are reported in Table4in terms of 𝑡 𝑐 .

Based on the results in Table

4

, the polynomial fittings of order 10 are selected. Polynomials of order 10 produce accurate and valid interpolations over a wide range of the parameter values. The validity interval 𝑡 𝑐 ∈ [0.06, 0.95] allows for solid volume fractions 1% ≤ 𝐵𝑉/𝑇𝑉 ≤ 99%.

Table 5 :

 5 Quality assessment of polynomial interpolation of the 𝐵𝑉 𝑇𝑉 ⁄ data.

	Polynomial order	R 2	RMSE	SSres
	2	0.9987964 0.017454 12.79248
	3	0.9999748 0.002527 0.26814
	4	0.9999973 0.000821 0.02829
	5	0.9999992 0.000442 0.00821

  fmincon is gradient-based, and it works on problems where the objective and constraint functions are both continuous and have continuous first derivatives. fmincon finds the minimum of a problem specified by 𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ] T and the objective function is 𝑓(𝑥) = ℛ 1 (𝑡 𝑐 , 𝑡 ℎ , 𝑡 𝑣 , 𝑡 𝑒 ), see equation (7).

				𝑐 𝑒𝑞 (𝑥) = 0
	min 𝑥	𝑓(𝑥) such that	𝑐 𝑖𝑛𝑒𝑞 (𝑥) ≤ 0 𝐴 𝑖𝑛𝑒𝑞 (𝑥) • 𝑥 ≤ 𝑏 𝑖𝑛𝑒𝑞	(11)
			{	𝑏 𝑙𝑜𝑤 ≤ 𝑥 ≤ 𝑢𝑝
	where x = [

Table 6

 6 reports the mean values of the residuals and the mean values and standard deviations of the relative errors of the geometric parameters and the symmetry classes.Besides, errors for the symmetry classes are relative to the corresponding average value of the 100 target microstructures. This approach avoids the occurrence of boundless and misleading large errors for the target microstructures with zero or nearly zero symmetry classes. Thus, the errors for the symmetry classes are defined as follows

	Errors for the geometric parameters are
	𝑒 𝑡 𝑐 = 𝑒 𝑡 𝑣 = 𝑡 𝑣 -𝑡 𝑣 𝑡 𝑐 -𝑡 𝑐 𝑡𝑎𝑟𝑔𝑒𝑡 𝑡 𝑐 𝑡𝑎𝑟𝑔𝑒𝑡 𝑜𝑡𝑎𝑟𝑔𝑒𝑡 𝑡 𝑣 𝑡𝑎𝑟𝑔𝑒𝑡	,	𝑒 𝑡 ℎ = and 𝑒 𝑡 𝑒 = 𝑡 ℎ -𝑡 ℎ 𝑡𝑎𝑟𝑔𝑒𝑡 𝑡𝑎𝑟𝑔𝑒𝑡 𝑡 ℎ 𝑡 𝑒 𝑡𝑎𝑟𝑔𝑒𝑡 𝑡 𝑒 -𝑡 𝑒 𝑡𝑎𝑟𝑔𝑒𝑡 ,	.		(17)
	𝑒 𝑐 𝑖𝑠𝑜 =	𝑐 𝑖𝑠𝑜 -𝑐 𝑖𝑠𝑜 𝑡𝑎𝑟𝑔𝑒𝑡 𝑐 𝑖𝑠𝑜 𝑡𝑎𝑟𝑔𝑒𝑡 ̅̅̅̅̅̅̅̅̅	,	𝑒 𝑐 ℎ𝑒𝑥 =	𝑡𝑎𝑟𝑔𝑒𝑡 𝑐 ℎ𝑒𝑥 -𝑐 ℎ𝑒𝑥 𝑐 ℎ𝑒𝑥 𝑡𝑎𝑟𝑔𝑒𝑡 ̅̅̅̅̅̅̅̅̅	,	(18)

Table 6 :

 6 Mean values of the residuals and errors for the geometric parameters and the symmetry classes.

  Table 6 shows that the SQP optimization produces the best results for the geometric parameters. SQP errors for the geometric parameters diminish with the relaxation of the 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 , maximum errors are around 1% for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 0 and they reduce to less than 0.1% for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 5%. Maximum errors are for 𝑡 𝑐 and 𝑡 ℎ . Standard deviations are, in general, within a few percent. Symmetry class errors behave almost independently of 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 . The result for 𝑐 𝑖𝑠𝑜 is very accurate, but in contrast, errors and standard deviations for 𝑐 ℎ𝑒𝑥 , 𝑐 𝑡𝑒𝑡 and 𝑐 𝑜𝑟𝑡 are very large. Mean values of 𝐵𝑉 𝑇𝑉 ⁄ for the optimal microstructures are almost coincident to the target values, with almost no effect of 𝐵𝑉 𝑇𝑉 ⁄

𝑡𝑜𝑙 ; mean values for the relative error

Table 7 .

 7 Errors are reported relative to the individual elasticity coefficients, It is observed that coefficients 𝐶 12 , 𝐶 13 and 𝐶 11 have the maximum individual relative errors, which range from 80% to 210%; however, when evaluated relative to ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖, maximum errors do not exceed 10%.

			𝑒 𝑖𝑗 𝐶 =	𝑡𝑎𝑟𝑔𝑒𝑡 𝐶 𝑖𝑗 -𝐶 𝑖𝑗 𝐶 𝑖𝑗 𝑡𝑎𝑟𝑔𝑒𝑡	,	(22)
	and relative to the elasticity matrix norm
			𝑒 𝑖𝑗 ‖𝐶‖ =	𝑡𝑎𝑟𝑔𝑒𝑡 𝐶 𝑖𝑗 -𝐶 𝑖𝑗 ‖ℂ 𝑡𝑎𝑟𝑔𝑒𝑡 ‖	.	(23)
	Error	Relative errors for elasticity matrix coefficients
	𝑒 𝑖𝑗 𝐶	𝐶 11 mean 0.80 -0.04 -0.24 0.34 0.86 2.15 -0.28 0.05 -0.21 𝐶 22 𝐶 33 𝐶 23 𝐶 13 𝐶 12 𝐶 44 𝐶 55 𝐶 66
		SD	1.06 0.34 0.26 0.66 0.83 1.49 0.38 0.54 0.81
	‖𝐶‖ 𝑒 𝑖𝑗	mean 0.10 -0.02 -0.11 0.01 0.04 0.08 -0.04 0.00 -0.01 SD 0.13 0.08 0.14 0.03 0.03 0.03 0.04 0.03 0.04

Table 7 :

 7 Relative errors for the elasticity matrix after the PS optimization.

Table 8 :

 8 Geometric parameters and 𝐵𝑉 𝑇𝑉 ⁄ of the mimetic bovine samples.

				parameters		𝐵𝑉 𝑇𝑉 ⁄
		𝑡 𝑐	𝑡 ℎ	𝑡 𝑣	𝑡 𝑒	[%]
	1	0.080 0.674 0.949 0.763	25
	2	0.187 0.229 0.931 0.600	38
	3	0.217 0.217 0.521 0.890	30
	4	0.080 0.391 0.718 0.727	20
	5	0.100 0.476 0.497 1.256	21
	Sample	𝑒 𝑐 𝑖𝑠𝑜 ,	Errors for symmetry classes 𝑒 𝑐 ℎ𝑒𝑥 , 𝑒 𝑐 𝑡𝑒𝑡 ,	𝑒 𝑐 𝑜𝑟𝑡 ,
	1	0.002	0.011	-0.076	0.018
	2	0.024	0.166	0.332	-0.696
	3	0.001	0.006	-0.049	0.013
	4	-0.002	-0.010	0.074	-0.019
	5	-0.001	-0.006	0.039	-0.004

Table 9 :

 9 Errors for the symmetry classes of the mimetic bovine samples.

	Sample			Relative errors for elasticity matrix coefficients
		𝐶 11	𝐶 22	𝐶 33	𝐶 23	𝐶 13	𝐶 12	𝐶 44	𝐶 55	𝐶 66
	1	0.034 0.016 -0.008 -0.039 -0.047 0.118 -0.028 -0.082 -0.032
	2	0.060 -0.016 -0.018 0.007 0.012 0.065 -0.007 0.015 -0.072
	3	-0.004 0.041 -0.035 0.005 0.003 0.048 -0.016 0.009 -0.054
	4	0.020 0.034 -0.048 0.029 0.014 0.087 -0.023 -0.001 -0.092
	5	0.167 -0.168 -0.039 -0.101 0.047 0.143 -0.062 0.003 -0.009

Table 10 :

 10 Relative errors for the elasticity matrix of the mimetic bovine samples.

ACKNOWLEDGEMENTS

This work has been supported by projects PIRSES-GA2009_246977 "Numerical Simulation in Technical Sciences" of the Marie Curie Actions FP7-PEOPLE-2009-IRSES of the European Union and by the PICS project "Modeling and Simulation in Multidisciplinary Engineering" MoSiMe funded in the framework of the CAFCI call by CONICET (Argentina) and CNRS (France). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 682001, project ERC Consolidator Grant 2015 BoneImplant).

CONFLICT OF INTEREST

The authors declare that they have no conflict of interest.

APPENDIX Human Samples

Erreur ! Source du renvoi introuvable. depicts the residual ℛ 1 and the symmetry class errors (18) of the SQP optimizations as functions of 𝐵𝑉/𝑇𝑉. Figure 11(a) shows that ℛ 1 and its dispersion diminish sharply with 𝐵𝑉/𝑇𝑉, from 0.5 ≲ ℛ 1 ≲ 3 for 𝐵𝑉/𝑇𝑉 < 7.5% to 0.1 ≲ ℛ 1 ≲ 0.6 for 𝐵𝑉/𝑇𝑉 > 32%. The relaxation of the 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 constrain does not result in significant improvements in ℛ 1 (results not reported here for analyses performed for 𝐵𝑉 𝑇𝑉 ⁄ 𝑡𝑜𝑙 = 10% exhibit the same behavior). Results in Erreur ! Source du renvoi introuvable.(b) allow to observe that 𝑐 𝑖𝑠𝑜 and 𝑐 ℎ𝑒𝑥 are systematically over and underestimated, respectively; their mean errors are 𝑒 𝑐 𝑖𝑠𝑜 ̅̅̅̅̅ = 0.20 and 𝑒 𝑐 ℎ𝑒𝑥 ̅̅̅̅̅̅ = -0.74. Errors tend to decrease with 𝐵𝑉/𝑇𝑉, the only exception is that of the tetragonal symmetry, which can attain values over the range -5 < 𝑒 𝑐 𝑡𝑒𝑡 < 12. Mean relative errors for the tetragonal and the orthorhombic symmetry classes are 𝑒 𝑐 𝑡𝑒𝑡 ̅̅̅̅̅̅ = 1.74 and 𝑒 𝑐 𝑜𝑟𝑡 ̅̅̅̅̅̅ = 0.18, respectively.