Mechanical bone growth stimulation by magnetic fibre networks obtained through a competent finite element technique

Fibre networks combined with a matrix material in their void phase make the design of novel and smart composite materials possible. Their application is of great interest in the field of advanced paper or as bioactive tissue engineering scaffolds. In the present study, we analyse the mechanical interaction between metallic fibre networks under magnetic actuation and a matrix material. Experimentally validated FE models are combined for that purpose in one joint simulation. High performance computing facilities are used. The resulting strain in the composite’s matrix is not uniform across the sample volume. Instead we show that boundary conditions and proximity to the fibre structure strongly influence the local strain magnitude. An analytical model of local strain magnitude is derived. The strain magnitude of 0.001 which is of particular interest for bone growth stimulation is achievable by this assembly. In light of these findings, the investigated composite structure is suitable for creating and for regulating contactless a stress field which is to be imposed on the matrix material. Topics for future research will be the advanced modelling of the biological components and the potential medical utilisation.

Scientific RepoRts | 7: 11109 | DOI: 10.1038/s41598-017-07731-6 also non-metallic 26,27 fibre networks have been tested for their ability to improve bone cement mechanics. The exhibited properties document very clearly improvements from the composite structure for this orthopaedic application. In non-static or smart applications of fibre-matrix composites, material properties are amended during the application process. Magnetic nanoparticles combined with cellulose fibres are in the discussion for the design of data storage applications and magnetographic printing or filtering 14,28 . This work has been extended also to incorporate bacterial cellulose 29,30 . In a different context, magnetic fibre networks could be applied to design heat exchangers of variable drag 31,32 . The assembly of this present study has been discussed for its suitability as scaffold in tissue engineering 33 . In this case, the purpose of the fibres would be to deliver under magnetic actuation a mechanical stimulus for the enhanced growth of bone cells. Various other fibre network materials are also investigated for their application as tissue scaffolds 29,34 . The delivery of a mechanical stimulus from scaffolds for influencing cells and their development is an on-going research topic on tissue level and on organ level [35][36][37] . For bone cells, it is known that cyclical loadings at 1 Hz of the strain magnitude 0.001 are beneficial for the tissue's remodelling behaviour 38 . A beneficial effect of a stand-alone magnetic field on bone growth is equally known 39,40 .
Yet, the mechanical interaction between fibre phase and matrix phase in fibre-matrix composites is so far only insufficiently understood. The purpose of the present study is to make a contribution to this field and to investigate the mechanical response of a matrix inserted in the void phase of a fibre network under magnetic actuation. Previous work has already analysed for this specific material simplified single-fibre geometries under magnetic actuation 15 . Or global values have been predicted analytically for fibre assemblies 41,42 . This present study investigates complete fibre network geometries and analyses the matrix strain on local level. The results of this study are obtained by means of experimentally validated finite element (FE) models and lead to an analytical model for the local matrix strain in cube shaped fibre networks under magnetic actuation. The presented results were part of a dissertation project at the University of Cambridge 19 .

Mathematical notation
The following notation is used in this document. Scalars are given as x, vectors as x , and 2 nd rank tensors as x . The vector product "×" and dot product "·" are applied. Relations which are greater-than and approximately equal are indicated as "≳". All fibre network nodes n of one sample i form the node set N i . The total number of nodes in that sample i is defined as N i . Equality between two volumes is written " ⇔ ", a volume containing a volume subset is given as "⊂". The operator for geometry definition by the intersection of two bodies is "∩". 3D spheres are specified for midpoint and radius as O (midpoint,R) . Their respective boundary is written o (midpoint,R) = ∂O (midpoint,R) , geometric cube faces are written F.

Modelling Methods and Material
For the investigations of the present study, a FE model (Fig. 1) of linear elasticity is developed. It describes the matrix strain under magnetic actuation in interaction with the fibre material. The quantitative evaluation of local fibre network density and local strain fields follows the scheme of Eqs 10 to 13 and leads to the analytical model presented by this study (Eq. 16).
Material samples and fibre network geometry extraction. FE modelling of the fibre network follows steps of an approach based on beam theory 12 . Three cube shaped fibre network samples are used in this present study ( Fig. 1a and Table 1). This allows the investigation of three different values of fibre volume fraction f: 10, 15, and 20%. The influence of f on the mechanical response is of particular interest.
These networks are produced by N.V. Bekaert S.A. (Belgium) from American Iron And Steel Institute (AISI) 316 L (d ≈ 40 μm) or by Nikko Techno Ltd. (Japan) from AISI 444 steel fibres, assembled as fibre stacks, compressed, and sintered to fibre mats. During the sintering step, inter-fibre bonds form. The sample sections are cut by electronic discharge machining from the sintered mats and computed tomography (CT) scans are acquired by General Electric (Germany) for a resolution of res = 7.75 μm. Due to the available computing capacities, the sample cube volume is set to V = (775.00 μm) 3 . In the CT scans, this sample cube side length is the equivalent of 100 pixels.
A skeletonisation algorithm [43][44][45] is applied and returns as output the medial axis paths of the fibre bodies contained in the CT data. Resulting from the chosen CT scan resolution res = 7.75 μm, the location of all medial axes is defined inside the sample cubes on a [7.75 μm · 100] 3 -grid. Nodes defining network medial axis location are referred to as network nodes in the following. Their total number per sample N i is known to increase for greater f 17,18 and is given in Table 1 for the three network samples.

FE meshing.
While previous studies have simulated fibre networks for a void inter-fibre space 9,12 , the present study combines a fibre phase V f and a matrix phase V M in the sample volume V. The meshing of both phases and the defined interaction rules are discussed in the following (Fig. 1c). The mesh is implemented for the FE solver Abaqus 6.13 46 (Tables 2 and 3). By their definition in the present study, V M equals V while V f is only a subset: Eq. 1 simplifies the volume description by disregarding the marginal volume of fibres whose medial axis is located at a distance to the surface S of less than d/2.
Matrix phase. A mesh of regular hexahadra elements (C3D8R) fills V. These hexahedra represent V M . For simplifying the interaction with the fibre bodies, the hexahedra side length L hex is set to the CT scan resolution res = 7.75 μm = 1 pixel. The number of hexahedra in V follows as N hex = 100 3 , their volume as V hex = (7.75 μm) 3 .
The values of matrix stiffness E M are chosen as it would be applicable for three stages in bone tissue development: collagenous bone 47 , granulation tissue 48 , and immature bone 49 . The value for immature bone was chosen from the early phase of bone development. The value for collagenous bone is slightly greater than found in the literature for enhancing the FE solver convergence.
Beam assembly. The fibre medial axes are transferred in the present study into beam assemblies. Two Timoshenko beams 50-52 are implemented. Results for linear interpolation (B31) and quadratic interpolation (B32) are analysed and compared to predictions from previous work 12     of ν = 0.3. The sintered inter-fibre joints are implemented by torsional springs (CONN3D2). Their stiffness is simulated as K joint = sE f A f with s = 5 μm. This value of s has been found to a pproximate results from experimental measurements 12,17 . The number of network nodes in Table 1 refers to exclusively those nodes which define beam start points and end points. The implementation for the FE solver requires additional technical nodes. Those purely technical nodes provide additional integration points for B32 elements, indicate the normal relative to the beam axis, and one is needed for each CONN3D2 element.
Interaction rules. The implementation of the interaction between V fibre and V matrix (Fig. 1c) is based on a previously published model for osteosynthesis applications 53 . The regular hexahedra and the beam mesh share the 100 3 -node grid which is defined by the CT scan pixels. Three main simplifications are made in the modelling: • V f and V M overlap, a double section assignment exists for up to max(f) = 20% of V.
• The interaction between both phases is reduced to the shared node grid.
• The mechanical interaction is set to identical displacement u at each node.
This means that more complex mechanical interaction as it appears in biological tissue between cells and scaffold is not considered at this point. These simplifications are made while, as in the original model, the following holds:

FE simulation and BC.
In the case of a void inter-fibre space, the Cauchy stress tensor σ is defined for a fibre network consisting of V f and the void matrix phase V void In this study, the relationship is modified for V M . Eq. 4 defines now the stress fields in the fibre-composite as: This relationship might seem obvious. It is mentioned explicitly at this point because Eq. 4 is required for the evaluation of imposed matrix strain with regard to the magnitude of E M and with regard to the respective magnetic actuation.
Magnetic actuation and BC. In the simulations, a magnetic induction vector B is imposed (Fig. 1b). The ferromagnetic mechanical response of the fibres is modelled as moment vector τ B . It is calculated for each beam element of length L f under the assumption of complete magnetisation parallel to the beam axis 15 :  In Eq. 9, τ M of all fibres in the respective sample is added into the equilibrium condition of moments. The actuation model is applied for M s = 1.6 MA/m. This value of M s has been experimentally validated for ferritic AISI 444 15 . The geometries of the network samples were obtained for austenitic AISI 316L. As part of this study's modelling assumptions, the network is simulated for the hypothetical M s of AISI 444. It has been shown that the geometry of the AISI 316L networks can similarly also be produced by AISI 444 fibres 17,18 . The imposed magnetic is orientated along the x-axis. At the time of this study's design, experiments were carried out at the University of Cambridge in which this type of in-plane actuation was tested.
The sample surface S of the sample volume V is defined by six quadratic cube faces, two of them perpendicular each to one of the three geometric axes. The cube face F x , marked in red in Fig. 1a, is used for the definition of the boundary conditions (BC). All three translational degrees of freedom (DOF) of beam elements located on the cube face F x are kinematically constrained up to a depth of h BC = 77.50 μm into the material. This value of h BC was applied as it provides a match to the experimental stiffness magnitude 12,17 .
Equilibrium conditions. The FE solver is run under the assumption of linear elasticity and for the equilibrium conditions of forces (Eq. 6) and moments (Eq. 8). It is required that σ = div( ) 0 is fulfilled. Body force per volume f is set to 0 in the present study; gravitational forces being a typical example for f . Together with the surface traction vector t, it defines the equilibrium of forces. Surface traction itself is defined by σ and the outer-normal vector: As the magnetic response is modelled as moment vector τ B , the general condition for equilibrium of forces can be simplified for the present study. Considering that = f 0 it follows that: Relative to the origin of point vector x , the general form of the equilibrium of moments 55 is defined: The simplification of Eq. 6 by = f 0 can be repeated for Eq. 8. The sum of τ B which is imposed on the fibre volume V f is added. The modified equation representing the present study and considering the sum of imposed τ B follows as: The Lagrangian reference frame 56 is used for Eqs 6 to 9. The alternative Eulerian reference frame can be advantageous for the modelling of fluids 57,58 . During the study, the solver runs into up to eight singularities and one zero-pivot (Table 1). For the facilitation of the solving step, singularities and zero-pivots are constrained for their respective DOF. Due to the size of the beam mesh (>10 3 ) this number of additionally constrained nodes is negligible with regard to the presented mechanical analyses.
Network node counting and strain field averaging methods. Eqs 10 and 11 are applied in this present study for further analyses of local fibre network density. For that purpose, a centre point μm is defined on the constrained cube face F x (Fig. 2a). The function y (R) stands for the relative share of nodes in sample volume V which is located on the sphere boundary o (CP,R) of radius R.
, 15%, 20%}: samples of 10, 15, and 20% (11) Y (R) stands for the integral of y (R) over R and counts network nodes of V inside the sphere O (CP,R) . y (R) and Y (R) lead to g (R) which is the first of two applied strain averaging methods. Both strain averaging methods, g (R) and f (R) , quantify the local strain obtained in the matrix phase V M . Two equivalent strain measures, von-Mises strain measure ε v.M. and maximum principal strain ε Pmax 59,60 , are calculated from the 3D strain fields of is generally recommended for ductile material, ε Pmax is of greater precision for brittle materials 61 . The material behaviour type of bone has been shown to depend on the deformation rate dε/dt 62 . g (R) quantifies the matrix strain on the boundary o (CP,R) where it intersects with V ( Fig. 2a and Eq. 12). f (R) places in a first step a sphere O (n,R) on each matrix node n and averages the matrix strain within O (n,R) . In the second step, the average is obtained across the entire node set N i (Fig. 2b and Eq. 13). g (R) quantifies the matrix strain depending on the distance to the constrained cube face. f (R) quantifies the matrix strain depending on the average distance to the fibre network structure.
. . . Hardware. For the solving step, the FE code is run on the Cambridge High Performance Computing Cluster Darwin. The cluster provides a total of 9600 cores by 600 quad server Dell C6220 chassis. Two 2.60 GHz eight core Intel Sandy Bridge E5-2670 processors with sixteen cores in total form one node with 64 GB of RAM (4 GB per core) 63 .

Results
The matrix phase of the meshed fibre-matrix composite (Fig. 1) is analysed for its deformation. Results about deformation patterns and deformation magnitude under magnetic actuation are presented. Distinctive patterns in the local matrix strain distribution are discovered. The results lead to an analytical model for the matrix strain magnitude in cube shaped fibre-matrix composites. A particular focus is the influence of the local fibre network density.
Local fibre network density. As explained above, the analysed fibre material is a statistical material. This requires for mechanical modelling knowledge about the distribution of fibre density and about the distribution of fibre orientation. In particular, this knowledge is of importance for the discussion of the local deformation patterns discovered in the present study. Previous analyses have shown that the main fibre orientation direction lies in the xy-plane; in the xy-plane itself, no preferred direction of orientation can be seen 17 . Figure 3 and Table 1 of the present study analyse now further the local distribution of material density across the three samples by quantification of the network node distribution. The density of network nodes relates directly to the fibre density in V.
The main result is that the variation of node density along the three Cartesian axes is far greater than the density variation relative to the sphere O (CP,R) . The network nodes are counted per 2D pixel slice of a 25 3 -grid along each of the three Cartesian axes. Table 1 shows the obtained numerical average, the median, and the standard deviation. A visible deviation exists in each of the three samples along each Cartesian axis. The maximum standard deviation relative to the average exists in min(f) = 10% along y and z. These two findings are very interesting as their cause lies in the compression step during manufacturing. Future network design might be able to exploit this further. However, this would make necessary a degree of control about production parameters greater than the one available in today's process. The corresponding distributions are plotted in Fig. 3a to c together with each sample's average. These three plots do not exhibit distinct local sample sections of outstanding greater/lesser material density. However, the extent of density variation along the three Cartesian axes becomes visible in these plots. The plot of Sample-10% shows clearly the greatest density variation relative to the sample average. Since the magnetic induction is imposed on the matrix material by the fibre network (Eq. 5), the rather uneven material density distribution along the three axes leads to local variations of obtained matrix strain (Figs 4 and 5b).
The share of nodes located on/in O (CP,R) is plotted in Fig. 3d and e as function y (R) and Y (R) (Eqs 10 and 11). Both, y (R) and Y (R) , exhibit a variation of the density distribution along R which is of far lesser magnitude than the one along the three Cartesian axes. This holds for each of the three samples. In Fig. 3d, y (R) is plotted. The three samples follow the same pattern of three distinct sections along R which can be seen in the distribution: • R ≤ 387.5 μm: y (R) increases exponentially. O (CP,R) corresponds to the smallest sphere drawn in Fig. 2a and is entirely contained in V. • 387.5 μm < R ≤ 775 μm: O (CP,R) exceeds the extensions of V at the cube's sides and y (R) decreases, in a good approximation linearly. • 775 μm < R: O (CP,R) intersects with V only at the four cube corners opposite to cube face F x as drawn for the greatest sphere in Fig. 2a. y (R) drops sharply to 0.
In Fig. 3e, it can be seen that the value of Y (R) increases for greater R continuously at different rates. A node share of 100% is reached when the entire V is contained in O (CP,R) , i.e R > 950 μm. Very important for the context of this study, the patterns are independent of f (i.e. it is obtained in each of the three samples). The findings of Fig. 3d and e are new for this specific material and are used in this present study for deriving the analytical model of the matrix strain magnitude (Fig. 6b and Eq. 16).

Matrix strain magnitude and distribution.
The FE assembly is simulated for the three sample geometries. The obtained magnitude of matrix strain and its distribution inside the sample, as well as, its statistical distribution are discussed now in this section. Deformation plot. Figure 4 contains the deformation plots of Sample-15% under magnetic actuation for all three matrix materials defined in Table 3. The green range on the colour bar is scaled to the magnitude 0.001 because of the particular interest with regard to applications including biomechanical bone growth stimulation 38 . Each of the nine plots contains locally this specific magnitude of matrix strain. The required magnitude of B is within the experimentally achievable range. Up to B = 1.00 can be seen as practical. This finding is of importance because it suggests the general suitability of this specific network material as scaffold for bone growth stimulation under magnetic actuation. Local maxima and minima are visible on S. For an increase from B = 0.25 to 1.00, the extension of the green areas increases, spreading further into the material. Greater E M reduces them respectively. The quantitative analysis of the matrix strain magnitude and distribution follows in Fig. 5.
Strain magnitude. Figure 5a plots the strain magnitude when averaged across the sample volume V for the previously shown interesting value of B = 1.00. The mechanical response is scaled linearly by B due to the modelling assumption of linear elasticity. The values of ε Pmax for brittle behaviour are marginally smaller than the values of ε v.M. for each value of f, and for each E M . This relationship is also predicted in the literature 59 and holds independently of the Timoshenko beam (Eq. 15). Lesser E M increases the matrix strain non-linearly. Of interest is the result that the average matrix strain clearly exceeds 0.001 in the case of collagenous bone and granulation tissue. At this early stage of bone tissue development, lesser values of B are sufficient. For immature bone and B = 1.00, average matrix strain between 10 −4 and 0.001 is obtained. A greater density of fibre network material increases the matrix strain for the given size of V. It is known for this network material that its mechanical response is influenced at this scale by a pronounced size effect 12 . Whether and how this size effect also influences the matrix strain under magnetic actuation is at this point one very worthwhile topic of future research.
The two Timoshenko beams lead to almost identical results for the matrix strain magnitude (Table 4). B32 is known to marginally reduce the network's mechanical stiffness 12 :  The obtained difference Δ is below the range of the magnitude of Δ is nearly independent of B. The obtained variation for B is in the range of the computing precision. The influence of f is not definite. For practical considerations in future prototype studies, the difference between both Timoshenko beams is negligible, especially with regard to other errors such as the one of the skeletonisation step. It is interesting that the difference increases for greater E M . One possible reason is that the stiffness requirement between matrix and fibres (Eq. 2) is better fulfilled for lesser E M . For numerical considerations, a further investigation of this model aspect through additional sample geometries will be highly worthwhile.
Spatial distribution. The spatial distribution of matrix strain (plotted on a 5 2 -grid and averaged along the z-axis) shows a comparable pattern between the three samples including at the same time local variations (Fig. 5b). The main finding is that the matrix strain magnitude is not uniform across the samples under magnetic actuation. For uniaxial mechanical actuation, the samples show a very uniform deformation 12 . Under magnetic actuation in this present study, the minimum matrix strain can be found in the proximity of the kinematically constrained cube face F x at min(x). Matrix strain magnitude increases with proximity to the non-constrained cube faces. In each of the samples, a maximum strain peak is located at the grid cell of max(x) and max(y). Inside the sample V, maxima of matrix strain exist too in each case. In the case of Sample-15%, a region of greatest strain on the cube's top face visible in the plots of Fig. 4 corresponds to the local maximum in Fig. 5b. Consideration of the results in Fig. 5b will be mandatory for future prototype designs. The strain magnitude and its location clearly depend on the applied BC and local variations exist. The design of BC imposed on the assembly will have to be considered in future studies. This is found to be a novel and reproducible feature of this specific material.    Matrix strain relative to fibre network and CP. The by far most relevant findings for the intended application purpose of this present study are shown in Fig. 6. Based on the strain distribution functions g (R) and f (R) (Eqs 12 and 13), it is possible to quantify and to analytically localize areas of greater or lesser matrix strain magnitude inside the material. Figure 6a plots f (R) which is the average matrix strain in O (n,R) depending on R averaged for all network nodes n of the sample. The obtained results show that R has a strong influence on f (R) . A specific pattern is discovered independently of f, E M , or B. In the close proximity to the fibre structure (R → 0 + ), f (R) exhibits in each sample one global peak which exceeds the sample average by several magnitudes. For increasing R ∈ [200 μm, 800 μm], f (R) drops below the sample average. In this range, f (R) is under the influence of the low-strain regions inside the samples and along F x (Fig. 5b). For values of R > 800 μm, the regions of greater matrix strain near the non-constrained cube faces influence f (R) stronger and it tends towards the sample average. This finding is particularly interesting for future research. One topic will be the mechanical interaction between fibres and the matrix at the interphase of these two. Biomechanical materials as used in the potential application of tissue engineering are characterized by non-linear mechanical response and additional parameters such as limited adhesion forces 64 . Those effects will be particularly important at the interface. Whether and how the predicted strain field changes remains to be seen.
The plots of Fig. 6b show g (R) , the matrix strain obtained on a sphere o (n,R) according to the averaging method Eq. 12. The obtained matrix strain increases for greater R and a local maximum is located at around R = 600 μm. This pattern is obtained for each matrix material in each sample. The individual differences of network geometry and material density in the samples do not change this overall pattern. The absolute maximum of strain is obtained in the non-constrained sample corners at R ≈ 900 μm for Sample-10% and Sample-20%. In the case of Sample-15%, the sample maximum is reached at R ≈ 850 μm. In each of the three samples, the obtained matrix strain drops afterwards. One difference between the samples is that this drop is much more pronounced in the case of Sample-15%. Low magnitude of matrix strain in the non-constrained corners of this sample is visible in the deformation plots (Fig. 4). It is worth mentioning that at this value of R only a fraction of the total 10 6 matrix hexahedra are located. The intersection of o (n,R) with V for R > 850 μm in the sample corners affects only ≈2% of all matrix hexahedra. This causes the obtained differences which are statistically not representative.
Based on the results of Fig. 6b, it is possible to derive for the first time a simple analytical model for the prediction of local matrix strain inside the sample's V. Eq. 16 approximates the obtained values by a power function: For R > 860 μm, the fit between Eq. 16 and the data declines sharply. This is why the regression analysis is run for R ≤ 860.25 μm. This means that 9,724 of 10 6 hexahedra are not included in Eq. 16, ≈1%. The matrix strain of the other 99% is modelled by Eq. 16. Table 5    specific pattern around CP and that it can be described by this simple analytical model. Differences can be seen between the three samples. The overall pattern repeats.

Discussion
We have presented above the results for the deformation of a fibre-matrix composite under magnetic actuation. The local matrix strain has been studied and quantified in depth for the first time. The results show for the investigated assembly distinct deformation patterns under magnetic actuation and conclusions about the influence of model parameters can be drawn. With regard to the general material response, it must be taken from now on into consideration that the matrix deformation is not uniform across the samples under magnetic actuation. This is in contrast to previous work which has been able to show that this particular network material exhibits under uniaxial mechanical actuation a very uniform deformation 12 .
One additional conclusion from this study is that BC have a direct influence on the resulting deformation field under magnetic actuation. For the cube shaped samples in this study, the local matrix strain reaches its minimum in the proximity of the constrained cube face. Local maxima exist inside the sample volume V. Maxima are located in the free corners of the sample cubes. In future studies, the applied BC will have to be explicitly considered. One major material parameter is the fibre volume fraction f. It can be concluded for this model parameter that the magnitude of matrix strain increases for greater f. Whether this trend will also be obtained for other sample volume V is one very interesting topic for future research. A pronounced size effect of the mechanical response is known for this network material 12 .
The statistical distribution contains in each sample one single global peak which shifts depending on matrix stiffness and magnetic induction. The conclusion is that controlling the strain magnitude of this incidence peak is mandatory for controlling the average and overall response in possible future applications. At the same time, also local sections of greater or lesser strain magnitude exist inside the matrix. The results for two derived distribution functions (Eqs 12 and 13) allow to gain novel insights about the localization of matrix strain magnitude inside the material. First, it can be shown that the matrix strain increases in each sample in the proximity of the fibre structure. Second, the distance to the kinematically constrained cube face has a strong influence on matrix strain magnitude and can be used for deriving a simple analytical model of matrix strain magnitude (Eq. 16). This model is valid for 99% of the matrix volume. The achieved fit for the regression analysis is a very good match with r 2 = 0.82 or greater. A remaining challenge are the free, unconstrained corners of the cube.
With regard to the purpose of this study to investigate the suitability of the fibre-matrix composite for possible applications as bioactive scaffold in tissue engineering, one very important result is obtained. From experimental work, it has been known since 1985 that a cyclical strain of 0.001 of the frequency 1 Hz is beneficial for the stimulation of bone growth 38 . This strain magnitude must be considered as design target for the scaffold application. The results in this study are showing that this magnitude is achievable by this investigated fibre-network composite assembly. In conclusion, it is possible to say that the principal suitability of this material for bone growth stimulation under magnetic actuation is confirmed. Future work will have to relax the simplifications made for the matrix model. For the application as tissue engineering scaffold, understanding of local strain at the interface between cells and steel is mandatory. The advanced version of this current model could integrate biological aspects such as the detaching of focal adhesion points build by bone cells on the fibre structure 64 . Also a dynamic matrix model including the active migration of bone cells through the fibre network 65 after the cell seeding would be a worthwhile aspect and of great interest for understanding bone development. The current model still relies on a static mesh.
In future work, one further effect must not be neglected in the evaluation of bone growth. It is known that magnetic fields alone already can stimulate bone growth. It is applied in the treatment of avascular necrosis or pseudarthrosis 39 . Designs for hip joint prostheses make use of that effect 40 . If benefits are achieved by the proposed assembly it will have to be made clear what the actual cause of the improvement is. The strain field imposed by the scaffold or the direct reaction of the tissue to a simple magnetic field? In the ideal case, it will be possible to raise synergies from both effects.