Computational Study on Interfacial Interactions between Polymethyl Methacrylate-Based Bone Cement and Hydroxyapatite in Nanoscale

: Polymethyl methacrylate (PMMA)-based bone cement (BC) is a key material in joint replacement surgery that transfers external forces from the implant to the bone while allowing their robust binding. To quantitatively evaluate the effect of polymerization on the thermomechanical properties of the BC and on the interaction characteristics with the bone ceramic hydroxyapatite (HAp), molecular dynamics simulations were performed. The mechanical stiffness of the BC material under external loading increased gradually with the crosslinking reaction occurrence, indicating increasing load transfer between the constituent molecules. In addition, as the individual Methyl Methacrylate (MMA) segments were interconnected in the system, the freedom of the molecular network was largely suppressed, resulting in more thermally stable structures. Furthermore, the pull-out tests using HAp/BC bilayer models under different constraints (BC at 40% and 85%) revealed the cohesive characteristics of the BC with the bone scaffold in molecular detail. The stiffness and the fracture energy increased by 32% and 98%, respectively, with the crosslink density increasing.


Introduction
Bone cement (BC) is a medical material used in joint replacement surgeries as a stabilizing anchor between the bone and the implant. Among the BC materials, poly (methyl methacrylate) (PMMA)-based organic BC is exceptionally light and flexible compared with inorganic cement materials, and its simple fabrication enables wide and universal applications. PMMA-based BC is largely composed of two parts: an MMA liquid and a PMMA powder [1][2][3]. In an actual surgical environment, these two parts are mixed in a bowl by stirring with a sterilized stick until a high-viscosity liquid is obtained. The resulting BC is applied to the surgical site using a syringe or hand, allowing the polymeric melt to infiltrate between the bone and the implant and to interdigitate into the bone network. From the mechanical and structural perspectives, BC acts as an intermediate medium that transmits various everyday external loads from the bone to the implant. In addition, because it allows a uniform distribution of the applied load on the entire surface, the disadvantages of metallic implants, including periprosthetic fracture or tissue damage, can be effectively prevented [4].
The binding characteristics at the interface of bone and BC were an important design parameter early on. Tozzi, G. et al. (2012) used finite element method (FEM) to derive the mechanisms of local damage and failure at the interface between bone and BC under compressive loading. The applied force was highly concentrated on the contact or partially 2 of 13 interdigitated region between bone and BC, which indicates that the yield elements of both materials were dominated in the corresponding field [5]. Simultaneously, in the fully interdigitated region, there was little yielding of both elements. In the study of Landgraf, R. et al. (2015), they considered the coupled simulation of computational fluid dynamics and FEM to describe the injection of BC into the cancellous bone structure and to investigate the impact of curling. The strategy allowed quantification of localized residual stress within the bone tissue due to the shrinkage of BC material as polymerization progressed [6]. While many studies including the abovementioned have emphasized the importance of the binding characteristics between bone and BC for the reliability of the resulting product in arthroplasty, it is rare to explore the origin of their interaction characteristics at the molecular level. Also, a quantitative estimation of bone-BC binding at the atomic level is especially important because the nano-and micro-porous structure of the bone scaffold with an enormous surface-to-volume ratio has a significant effect on the performance of the entire system.
Meanwhile, the molecular conformation and properties of the PMMA-based BC were determined by polymerization via free radical reactions between the constituent molecules. In many studies conducted to date, the effect of crosslinking on the thermomechanical properties of polymer materials (such as stiffness, thermal expansion, and glass transition temperature) was very apparent [7][8][9][10][11][12][13]. However, because many of these quantitative results have been reported for some specific thermosetting resins used in various industries (such as epoxies and polyimides), characterization of PMMA-based BC remains relatively poor. Moreover, for the BC, it is difficult to monitor quantitative changes in its physical properties according to the crosslinking rate variations, because such variations are dominated by atomic-level compositional details, making it difficult to identify in experimental approaches.
Another issue associated with polymer crosslinking is the nature of the interactions with other materials. Although the interaction behavior occurs in a very localized region of several nanometers, these interfacial interactions affect the response of entire porous structures in the skeletal system. In particular, the constituent interfacial molecules exhibit different structural and mechanical responses than those in the bulk regions, owing to the non-bonded interactions between materials. Therefore, interfacial effects are strongly coupled to variables that affect morphological features, such as crosslinking density [12,14], and a comprehensive design strategy is required.
The present MD work is the first to provide a comprehensive understanding of the interaction properties between bone and BC from the mechanical perspective and has two main objectives. The first objective was predicting the thermomechanical properties of the PMMA-based BC according to the degree of crosslinking. For this, BC microstructures under four different crosslinking conditions were prepared, and both a tensile-loading simulation and a cooling-down simulation were conducted. The second objective was characterizing the interfacial cohesive properties with the bone scaffold. The bilayer model of the BC and hydroxyapatite (HAp) that a key bioceramic material of the bone scaffold was considered, and the fracture toughness as well as the cohesive strength at the BC/HAp interface were quantitatively characterized throughout the pull-out tests.

BC Microstructures
Amorphous BC microstructures with various crosslinking ratios were prepared using all-atom MD simulations. The molecular structures used in this study are shown in Figure 1a. To describe the two components of the PMMA-based BCs, an MMA composed of four monomers and a PMMA chain in which 50 monomers were linearly connected were prepared. Initially, both MMAs and PMMAs were randomly positioned, with the weight fraction of 1:2 in the unit cell. To describe a continuous bulk system, periodic boundary conditions were employed along all three axes. The unit cell was constructed with an initial density of 1.15 g/cc and the potential energy of the uncrosslinked BC model was then minimized using the conjugate gradient method (cutoff at 10 kcal/mol Å). Subsequently, an equilibration period using the Verlet algorithm composed of the normal volume temperature (NVT) dynamics at 300 K for 1 ns and the normal pressure temperature (NPT) dynamics at 300 K and 0.1 MPa for 1 ns was allowed, with the time step of 1 fs. All interatomic behaviors of constituent atoms were simulated using a polymer consistent force field (PCFF) [15].
tween them was measured. Then, covalent bonds between the reactive points were newly created starting from the shortest distance, until the total crosslinking degree in the unit cell reached approximately 80%. Regarding the crosslinking reaction, when the crosslinking ratio reaches 80%, almost all molecules inside the unit cell participate in polymerization as shown in Figure 1c. To observe the changes in the internal microstructure and interaction characteristics according to the crosslinking ratio, BC models under four different conditions (0%, 20%, 50%, and 80%) were considered. Each selected BC model was further equilibrated by the NVT dynamics at 300 K for 1 ns, and by the NPT dynamics at 300 K and 0.1 MPa for 2 ns. The details of the equilibrated BC microstructures are listed in Table 1. In all cases, the unit cell density after the NPT ensemble was in the 1.07-1.08 g/cc range, which was confirmed to converge under conditions similar to the typical values observed experimentally [16][17][18]. Note that these modeling procedures and crosslinking mechanisms for polymer microstructures used in this study can also be found in our previous works [19][20][21].  One key aspect of interest in the present study was the dependence of the physical properties on the crosslink density of the BC. An illustration of the crosslinking procedure and the corresponding BC microstructure, for the crosslinking ratio of 80%, are shown in Figure 1b. First, reactive points were assigned to each component, and the distance between them was measured. Then, covalent bonds between the reactive points were newly created starting from the shortest distance, until the total crosslinking degree in the unit cell reached approximately 80%. Regarding the crosslinking reaction, when the crosslinking ratio reaches 80%, almost all molecules inside the unit cell participate in polymerization as shown in Figure 1c. To observe the changes in the internal microstructure and interaction characteristics according to the crosslinking ratio, BC models under four different conditions (0%, 20%, 50%, and 80%) were considered. Each selected BC model was further equilibrated by the NVT dynamics at 300 K for 1 ns, and by the NPT dynamics at 300 K and 0.1 MPa for 2 ns. The details of the equilibrated BC microstructures are listed in Table 1. In all cases, the unit cell density after the NPT ensemble was in the 1.07-1.08 g/cc range, which was confirmed to converge under conditions similar to the typical values observed experimentally [16][17][18]. Note that these modeling procedures and crosslinking mechanisms for polymer microstructures used in this study can also be found in our previous works [19][20][21].

HAp Substrate and HAp/BC Bilayer Models
To evaluate the interaction between the BC and bone scaffolds with respect to the crosslinking degree, an HAp substrate model was considered, as shown in Figure 2a. We considered HAp (110) as the reference interaction plane, for estimating the possible maximal binding capability, for the highest planar atom density [22,23]. The HAp model with the (110) surface had rectangular lattice parameters of a = 6.879 Å, b = 16.323 Å, c = 9.922 Å and α = β = γ = 90 • [24]. The HAp substrate model was constructed by replicating the lattice with 7 × 3 × 4 in each direction. For the accurate description of the HAp substrate's behavior as well as its surface, the additional parameters produced by Lin, T. et al. [25] based on the PCFF were adopted for the bond, angle, and non-bond interaction terms. Periodic boundary conditions were used only for the x and y directions for the surface relaxation. The HAp model was equilibrated with the NVT dynamics at 300 K for 1 ns, and with the NPT dynamics at 300 K and 0.1 MPa for 2 ns, prior to merging with the BC model. During the series of dynamics, the HAp model well preserved their initial geometry including the surface region; the resulting HAp substrate had the dimensions of 48.01 × 48.98 × 38.7 Å 3 . The HAp/BC bilayer models were prepared by placing the BC models with different crosslinking ratios atop the HAp model, with the spacing of 1 nm. Interactions between HAp and BC occurred under the NVT dynamics at 300 K for 2 ns. As a result, the gap between the two dissimilar materials disappeared naturally, and both came into a full contact with each other, as shown in Figure 2b. Note that during the construction of the HAp/BC bilayer models, a sufficiently thick vacuum layer with a thickness of 7 nm was inserted at the top of the BC to completely remove interactions beyond the unit cell boundary along the z direction. The BC and HAp models were built using commercial software, Materials Studio ® (Dassault Systèmes), and the interactions between the BC and HAp models during MD runs were observed using the open-source code LAMMPS (Sandia Lab.).

Effect of the Crosslinking Reaction on the Thermo-Mechanical Properties
Uniaxial tensile loading simulations were conducted for various temperatures, ranging from 300 K to 390 K. The BC microstructure models were equilibrated at each temper-

Effect of the Crosslinking Reaction on the Thermo-Mechanical Properties
Uniaxial tensile loading simulations were conducted for various temperatures, ranging from 300 K to 390 K. The BC microstructure models were equilibrated at each temperature for 1 ns, using the NPT dynamics, and then were frozen immediately at 0.1 K for 0.5 ns, using the NVT dynamics. The tensile simulation near zero temperature was performed for each BC model until the overall strain reached 0.15% at a constant strain rate of 10 6 s −1 . The stress tensor (σ) of the BC unit cell was estimated from the virial theorem, as follows [26]: where V and N are the system volume and the number of constituent atoms, m i and v i are the atomic mass and the velocity of atom i, and r ij and F ij are the relative distance and the atomic force between atoms i and j, respectively. Because all of the components in σ were extracted at 0.1 K, where the kinetic term was eliminated without changing the direction of the traction vector during the deformation, the stress was equivalent to the Cauchy stress [27]. In addition, the stiffness tensor (C) of the BC unit cell model was calculated as follows: where ε represents the applied strain tensor. The tensile simulation was performed for each direction, and the stiffness results were averaged assuming an isotropic medium. The thermoelastic behavior characteristics of the models, for four different temperatures, are shown in Figure 3a. At room temperature, an increase in the BC stiffness was clearly observed as the crosslinking density increased (from 5.77 GPa to 6.27 GPa). This result was straightforward because the load transfer between the BC constituent molecules became more efficient and dominant as the individual MMA liquid became polymerized in the system. Notably, the magnitude of elastic stiffness at 330 K is quite similar to the one observed at 300 K. The reinforcement on the stiffness with the progress of crosslinking reaction is also valid. The similarity between these two temperature conditions agrees very well with the results observed in the experiment [28,29]. At temperatures above 360 K, the stiffness of the BC decreased considerably, and the deterioration of stiffness became more pronounced with increasing temperature. The sharp decrease of the elastic stiffness of BC material at temperatures above 360 K is in agreement with the experimental observations [28,29]. Furthermore, it is interesting to find that the influence of the crosslinking density varies greatly with increasing temperature. This implies that the mechanical properties of the BC material after being stabilized in the body are significantly different compared to when the BC material is placed in a high exothermic state by a crosslinking reaction in the surgical environment. Meanwhile, the structural stability of the BC material can also be evaluated from the change in the free volume between polymer chains. Figure 3b,c shows the change in free volume of the BC unit cell model at each stage of temperature. The free volume was calculated based on the Connolly surface, where the probe radius was set as 2 Å. As shown in the results, there was a rapid increase in free volume at temperature of 360 K and above, which causes a softening of the BC material [28].
To evaluate the thermal stability of the PMMA-based BCs at room temperature, mean squared displacement (MSD) curves were obtained at 300 K using the following equation: perimental observations [28,29]. Furthermore, it is interesting to find that the influence of the crosslinking density varies greatly with increasing temperature. This implies that the mechanical properties of the BC material after being stabilized in the body are significantly different compared to when the BC material is placed in a high exothermic state by a crosslinking reaction in the surgical environment. Meanwhile, the structural stability of the BC material can also be evaluated from the change in the free volume between polymer chains. Figure 3b,c shows the change in free volume of the BC unit cell model at each stage of temperature. The free volume was calculated based on the Connolly surface, where the probe radius was set as 2 Å . As shown in the results, there was a rapid increase in free volume at temperature of 360 K and above, which causes a softening of the BC material [28]. To evaluate the thermal stability of the PMMA-based BCs at room temperature, mean squared displacement (MSD) curves were obtained at 300 K using the following equation: where R i ⃗⃗⃗ (t) indicates the position vector of atom at time . Figure 4a shows the MSD curves obtained for the initial 200 ps in the NPT ensemble. The MSD of the uncrosslinked BC model was especially large over the entire simulation time, and the profiles gradually decreased with the crosslinking ratio. The results suggest that as the MMA liquid became covalently connected to the adjacent MMA and/or PMMA, the translational freedom of the molecular networks became highly suppressed. The linear coefficient of thermal expansion (LCTE) was also calculated for each BC model. Each BC model was equilibrated at 390 K and 1 atm using the NPT dynamics for 1 ns, and then was gradually cooled down at a constant rate of 0.02 K/ps. For each level of temperature, the change in the specific volume of the unit cell was monitored, and the LCTE was approximated using the following linear regression [30]: where is the LCTE, while is the volumetric coefficient. In addition, and are the volume and the specific volume, respectively, while 0 and 0 are the volume and the density of the unit cell at room temperature, respectively. Similar to the results for the MSD curves, the LCTE value of the fully crosslinked BC decreased by approximately 20%, compared with the uncrosslinked BC, as shown in Figure 4b. In particular, the LCTE value for BC-80 (80.5 ppm/K) agreed well with the experimentally reported value (approximately 98 ppm/K, for the highest molecular weight) [31]. The strong covalent bonds among the BC materials improved the thermal stability by limiting the free movement of the molecular chains. In conclusion, the change in the crosslinking density with curing time clearly affected both the structural configuration and the thermal stability of the microstructures. Such significant changes in the thermomechanical properties of BC materials are particularly important for dealing with interactions with the bone scaffold (HAp), as will be discussed in the next section.

Density Distributions in the Interfacial Regions of HAp/BC Bilayer Models
To observe the change in the interfacial behavior with respect to the crosslinking ratio, the density distribution along the z direction was measured for each HAp/BC model, and the results are shown in Figure 5. Points in the density profiles were obtained with The linear coefficient of thermal expansion (LCTE) was also calculated for each BC model. Each BC model was equilibrated at 390 K and 1 atm using the NPT dynamics for 1 ns, and then was gradually cooled down at a constant rate of 0.02 K/ps. For each level of temperature, the change in the specific volume of the unit cell was monitored, and the LCTE was approximated using the following linear regression [30]: where α is the LCTE, while γ is the volumetric coefficient. In addition, V and v are the volume and the specific volume, respectively, while V 0 and ρ 0 are the volume and the density of the unit cell at room temperature, respectively. Similar to the results for the MSD curves, the LCTE value of the fully crosslinked BC decreased by approximately 20%, compared with the uncrosslinked BC, as shown in Figure 4b. In particular, the LCTE value for BC-80 (80.5 ppm/K) agreed well with the experimentally reported value (approximately 98 ppm/K, for the highest molecular weight) [31]. The strong covalent bonds among the BC materials improved the thermal stability by limiting the free movement of the molecular chains. In conclusion, the change in the crosslinking density with curing time clearly affected both the structural configuration and the thermal stability of the microstructures. Such significant changes in the thermomechanical properties of BC materials are particularly important for dealing with interactions with the bone scaffold (HAp), as will be discussed in the next section.

Density Distributions in the Interfacial Regions of HAp/BC Bilayer Models
To observe the change in the interfacial behavior with respect to the crosslinking ratio, the density distribution along the z direction was measured for each HAp/BC model, and the results are shown in Figure 5. Points in the density profiles were obtained with spatial steps of 0.2 Å, and the zero point of the z-axis indicates the surface of the HAp substrate. The condensed BC phases were laid approximately 3 Å away from the HAp surface after relaxation. The non-bond interaction between the BCs and the HAp substrate plays a dominant role in the vicinity of the interface, but its effect became much weaker as the distance from the interface increased. One notable point is that the peak magnitudes of the density profiles increased with increasing the crosslinking ratio (from 1.25 g/cc to 1.44 g/cc). In terms of the interactions between polymer-heterogeneous material, these results are quite different from those of other studies using epoxy resins, which reported a decrease in the peak magnitude with increasing crosslinking density [12,14,32]. This difference can be attributed to the structure of the polymer network generated during the crosslinking reaction. For a three-dimensional network structure including epoxy thermosets, physical interlocking events induced by the crosslinking reaction occur widely and homogeneously among the polymer matrix [33], which work as spatial constraints that inhibit interfacial effects [12]. On the other hand, for linear polymer structures (i.e., PMMA used in the current study), crosslinking events occur in inhomogeneous media and cause a sparse entanglement within the microstructure. Consequently, the crosslinking reaction on linear BC chains only reduces the diffusivity by removing MMA segments, leading to stable and improved interfacial interactions with the HAp substrate.
Another interesting aspect is that the oscillation of the BC density profile became clearer as the crosslinking reaction progressed. This was attributed to the fact that linear polymer chains adsorbed on the heterogeneous material tend to generate regular lamellartype structures [34]. Thus, as the diffusivity decreased and a more stable interaction state was achieved, the oscillation became more pronounced. on linear BC chains only reduces the diffusivity by removing MMA segments, leading to stable and improved interfacial interactions with the HAp substrate.
Another interesting aspect is that the oscillation of the BC density profile became clearer as the crosslinking reaction progressed. This was attributed to the fact that linear polymer chains adsorbed on the heterogeneous material tend to generate regular lamellartype structures [34]. Thus, as the diffusivity decreased and a more stable interaction state was achieved, the oscillation became more pronounced.

Cohesive Characteristics of the HAp/BC Bilayer Model
To evaluate the dependence of the BC/HAp cohesion on the crosslinking density, two types of the pull-out tests were conducted in this study: (1) a test in the opening mode and (2) a test in the sliding mode. For the opening-mode test, the upper part of the BC sample was fixed and the HAp sample was moved in the -z direction at a speed of 0.5 m/s, using

Cohesive Characteristics of the HAp/BC Bilayer Model
To evaluate the dependence of the BC/HAp cohesion on the crosslinking density, two types of the pull-out tests were conducted in this study: (1) a test in the opening mode and (2) a test in the sliding mode. For the opening-mode test, the upper part of the BC sample was fixed and the HAp sample was moved in the -z direction at a speed of 0.5 m/s, using the NVT dynamics at room temperature. During the pull-out process, the traction stress was obtained by tracking the force applied to the HAp substrate and dividing it by the surface area. The sliding-mode test was performed by moving the HAp substrate in the x or y directions, under the same BC constraint and under the same substrate velocity condition, while the traction was derived by averaging the responses for the different directions.
First, the traction-separation response was calculated, while 40% of the top of the BC was fixed to derive the cohesive properties throughout the entire system. As shown in Figure 6a, when the separations were applied in the direction perpendicular to the surface, a linear increase in the traction stress was observed. In the corresponding elastic region, many uniform-size nanovoids were created in the polymer network, whereas the nanovoids did not grow [35]. Especially, the magnitude of the traction stress was proportional to the occurrence of crosslinking, indicating improved cohesive properties. In the subsequent damage region resulting from the continuous descent of HAp, these nanovoids gradually integrated over separation, decreasing the traction stress exponentially. Note that because only the partial region of the BCs was fixed, the current cohesive response is a comprehensive result that includes not only the BC/HAp interfacial characteristics, but also the mechanical response of the BC. To separately examine the contribution of the interfacial region, the same pull-out test was performed with the top 85% of the BC fixed, as shown in Figure 6b. Because the 85% constraint on the BC excludes the variation of the mechanical properties by the curing of the BC material, only the interfacial response can be identified under this condition. Still, stronger traction was clearly observed during the interface opening, for high crosslinking ratios. In conclusion, the crosslinking reaction of  Similar analysis was performed for the sliding mode, and the results are shown in Figure 7. The associated cohesive parameters are listed in Table 4 (for the entire system) and Table 5 (for the interfacial region). Note that periodic boundary conditions were applied to the x and y directions of the unit cell, for describing the case in which the HAp substrate was constantly and permanently engaged with the BC material. Therefore, the obtained responses in the sliding mode are only valid in the elastic region where traction In the traction-separation relationship, the δ 0 and δ f parameters represent the elastic distance and the failure distance, respectively. In addition, K e f f is the cohesive stiffness in the loading direction, and G c is the critical energy required for fracture, which can be calculated from the area under the curve. The characterized cohesive parameters are listed in Tables 2 and 3. In all of the considered cases, K e f f and G c significantly improved with a progress in the crosslinking ratio (K e f f improved by 32%; G c improved by 98% for the system separation. In addition, for the interface separation, K e f f improved by 14% while G c improved by 9%). In the case of the entire system separation, an increase in δ f was also observed because the effective length of the polymer network increased owing to the high crosslinking ratio.
The relationship between HAp and BC mainly originated from the weak non-bonded interactions rather than physical bonds such as covalent anchoring. For this reason, the critical problems within the body such as loosening of the resulting product due to the slippery interface can occur if the interaction strength is not sufficiently high. The present results suggest that as the polymerization of the BC material progresses, it has a higher cohesive stiffness and fracture toughness for both constraint conditions, which enables effective stress transfer under the various loading environments.  Similar analysis was performed for the sliding mode, and the results are shown in Figure 7. The associated cohesive parameters are listed in Table 4 (for the entire system) and Table 5 (for the interfacial region). Note that periodic boundary conditions were applied to the x and y directions of the unit cell, for describing the case in which the HAp substrate was constantly and permanently engaged with the BC material. Therefore, the obtained responses in the sliding mode are only valid in the elastic region where traction is continuously growing, rather than fully describing the behavior of typical tractionseparation. As shown by the results, the response in the sliding mode exhibited the same tendency as in the opening mode. The traction stress increased linearly in the elastic region and their magnitude also depends on the progress in the crosslinking reaction, for both constraint conditions. Accordingly, the K e f f values increased up to 40% and 69%, respectively. The distance of the elastic region, δ 0 , was also proportional to the crosslinking density for the separation of the entire system, indicating the enhanced fracture toughness of the system. is continuously growing, rather than fully describing the behavior of typical traction-separation. As shown by the results, the response in the sliding mode exhibited the same tendency as in the opening mode. The traction stress increased linearly in the elastic region and their magnitude also depends on the progress in the crosslinking reaction, for both constraint conditions. Accordingly, the values increased up to 40% and 69%, respectively. The distance of the elastic region, 0 , was also proportional to the crosslinking density for the separation of the entire system, indicating the enhanced fracture toughness of the system.

Conclusions
In this study, the thermomechanical properties of PMMA-based BC and the interaction features with a bioceramic material, HAp, were investigated using all-atom MD simulations. The two types of pull-out tests enabled to quantitatively evaluate the cohesive properties of the bilayer system with respect to the crosslinking density.
The results suggested that increasing in the crosslinking reaction improved the mechanical behavior of the BC material by facilitating the load transfer between the polymer chains. Moreover, the reduction in the chain mobility through the polymerization of the MMA segment ensured a stable interaction with HAp as well as enhanced thermal stability. Furthermore, in the opening and sliding modes of pull-out tests, it was verified that the strength and toughness between HAp and BC materials improved with the degree of crosslinking, and presented higher K e f f and G c . From the two different fixing conditions of 40% and 85% of the upper part of the BC, it turn out that the cohesive characteristics take advantage of both the variation in the mechanical response of the BC and in the interfacial behavior of the BC/HAp in terms of the crosslinking reaction.
The current MD study is the first step toward exploring the interaction characteristics between biopolymer cements and the skeletal system on the atomic level from a mechanical point of view. The results can be used in conventional continuum approaches, such as the finite element method (FEM) with a cohesive zone model (CZM) [36]. Hence, one can consider the multi-scale modeling regime that combines the classical fracture criterion with the critical energy (G c ) collected in this study. The local stress concentration in the body of a patient under mechanical shock and subsequent crack initiation can also be covered.