The effect of collagen fibril orientation on the biphasic mechanics of articular cartilage.

The highly inhomogeneous distribution of collagen fibrils may have important effects on the biphasic mechanics of articular cartilage. However, the effect of the inhomogeneity of collagen fibrils has mainly been investigated using simplified three-layered models, which may have underestimated the effect of collagen fibrils by neglecting their realistic orientation. The aim of this study was to investigate the effect of the realistic orientation of collagen fibrils on the biphasic mechanics of articular cartilage. Five biphasic material models, each of which included a different level of complexity of fibril reinforcement, were solved using two different finite element software packages (Abaqus and FEBio). Model 1 considered the realistic orientation of fibrils, which was derived from diffusion tensor magnetic resonance images. The simplified three-layered orientation was used for Model 2. Models 3-5 were three control models. The realistic collagen orientations obtained in this study were consistent with the literature. Results from the two finite element implementations were in agreement for each of the conditions modelled. The comparison between the control models confirmed some functions of collagen fibrils. The comparison between Models 1 and 2 showed that the widely-used three-layered inhomogeneous model can produce similar fluid load support to the model including the realistic fibril orientation; however, an accurate prediction of the other mechanical parameters requires the inclusion of the realistic orientation of collagen fibrils.


Introduction
The main function of articular cartilage (AC) is to serve as a bearing material for the synovial joints. Its unique load-bearing and lubrication properties are due, in part, to its composition which includes a solid phase of collagen fibrils enmeshed with proteoglycans and a fluid phase (Lu and Mow, 2008). The high tensile stiffness of collagen considerably increases the compressive stiffness of AC by providing resistance to lateral expansion and thereby pressurising the interstitial water . Fluid pressurisation is believed to be one important reason that AC exhibits a very low-friction coefficient (Krishnan et al., 2004).
The distribution of collagen fibrils in AC is highly inhomogeneous. It is believed to be depth-dependent: oriented to merge with the articular surface in the superficial zone (SZ), perpendicular to the boundary between bone and cartilage in the deep zone (DZ), and randomly oriented in the middle zone (MZ) (Clark, 1990;Kaab et al., 1998). Moreover, it also appears to be location-dependent (Clark, 1991). For example, in the centre of the tibial cartilage, the DZ is thick and orderly, whereas the SZ and MZ are thin (combined contribution E5%) and ill-defined. In the periphery, the SZ and MZ layers are thicker (50%) and well developed (Clark, 1991). Such an inhomogeneous distribution of collagen fibrils may have important effects on the mechanics of AC.
The effect of the inhomogeneous distribution of collagen fibrils has been mainly investigated using finite element (FE) models, mostly using simplified methods, in which cartilage is divided into three zones to simulate the SZ, MZ and DZ. The boundaries between the zones have been determined either by data from the literature (Dabiri and Li, 2013;Li et al., 2000) or the T 2 relaxation time of magnetic resonance imaging (MRI) (Julkunen et al., 2008;Rasanen et al., 2013). In these studies, the orientation of the collagen fibrils in the SZ, which appears to be slightly oblique to the articular surface (Clark, 1990;Kaab et al., 1998;Meder et al., 2006;Wu et al., 2008), was assumed to be perfectly parallel to the articular surface; while in the DZ, it was assumed to be perfectly perpendicular to the cartilage-bone interface. Moreover, although the effects of the inhomogeneous distribution of collagen fibrils can also be investigated using the formulation proposed by Federico et al. (Federico and Gasser, 2010;Federico and Grillo, 2012;Federico and Herzog, 2008c), which can model the mechanical response of the statistically distributed collagen fibres in arbitrary directions across the cartilage layer, the published applications of this promising formulation are still limited (Tomic et al., 2014;Wu et al., 2016). Therefore, how the above simplifications affect the prediction of biphasic cartilage mechanics is still not clear due to the lack of direct comparisons with models considering a realistic orientation of the collagen fibrils.
A more accurate representation of the collagen fibril orientation in the FE models than the simplified three-layered orientation distribution will be location-dependent and therefore more inhomogeneous. Since considering the effect of fibrils in computational models considerably increases the fluid pressure Soltz and Ateshian, 2000;Soulhat et al., 1999) and the inclusion of the depth-dependent inhomogeneity of cartilage further increases the interstitial fluid load support at the articular surface (Dabiri and Li, 2013;Krishnan et al., 2003), it is reasonable to expect that the location-dependent orientation of collagen fibrils may change the interstitial fluid load support further. Moreover, the inclusion of the tensile effect of fibrils and the depth-dependent inhomogeneity of cartilage in computational models also substantially affects the cartilage mechanics, i.e., the distribution of the stresses and strains (Dabiri and Li, 2013;Krishnan et al., 2003;Li et al., 2000;Soulhat et al., 1999). Therefore, a more inhomogeneous collagen fibril distribution may predict considerably different stress and strain states. This in turn will affect the prediction of mechanically induced damage, which has been linked to the maximum shear stress (Atkinson et al., 1998;Eberhardt et al., 1991), tensile stress (Kelly and O'Connor, 1996), principal strain (Wilson et al., 2006), and shear strain (Wilson et al., 2006).
The development of MRI techniques has enabled the quantification of the orientation of collagen fibrils of AC. For example, based on the assumption that the diffusion of water in the direction of fibres is faster than in the perpendicular directions, the orientation of collagen fibrils can be quantified by the principal eigenvector of the diffusion tensor, which can be measured by diffusion tensor MRI (DT-MRI) (de Visser et al., 2008b;Filidoro et al., 2005;Meder et al., 2006). Recently, Pierce et al. (2013) investigated the inhomogeneities of the collagen fibre distribution derived from the DT-MRI data and the material properties on the interstitial fluid pressure of AC. However, a comprehensive understanding of the effect of the realistic orientation of collagen fibrils on the cartilage mechanics is still not clear, especially with regard to what extent the widely-used three-layered model can represent the realistic situation under time-dependent conditions. Therefore, the aim of this study was to investigate the effect of the realistic orientation of collagen fibrils on the biphasic mechanics of cartilage using a model based on DT-MRI data.

DT-MRI data
A cylindrical cartilage-on-bone sample (6 mm diameter) was extracted from the patellofemoral groove of a slaughtered cow (aged approximately 18 months), placed in phosphate-buffered solution (PBS) and frozen prior to use. Measurements were performed after thawing in PBS at approximately 20°C. The cartilage thickness was approximately 1.2 mm.
Measurements were carried out on a Bruker Avance II, vertical bore spectrometer, operating at 9.4 T. The spectrometer was equipped with a 25 mm diameter RF coil and actively shielded gradients which were supplemented with pre-emphasis. With the imaging parameters as given below, no eddy current effects were observed in the images. A standard, spin-echo based, single-echo, diffusion-tensor MRI sequence (Bruker, ParaVision 4) was used with full k-space sampling. Ten unique diffusion gradient directions were used according to an icosahedral scheme equivalent to the ICOSA10 design of Hasan et al. (2001) (see Appendix A for details).
Sequence timing parameters were T E ¼55.9 ms, T R ¼ 2000 ms. A single slice was acquired with a thickness of 2 mm. A matrix size of 256 Â 256 was used, giving an isotropic in-plane pixel size of 70 mm. For each pixel in the images, diffusion tensor components were calculated using an unweighted multilinear regression (using in-house MATLAB code) against b-values calculated from the diffusion gradients and interactions with the spoilers and slice-selection gradients (negligible contribution). In addition, linear background gradients were also fitted by taking into account the interaction between such gradients and the diffusion and spoiler gradients. Thus, ten parameters were fitted within each pixel: a constant term, six diffusion components, and three vector components. More details of the calculation of the diffusion tensor can be found in Appendix B.
Only the tensor components were used in this work and, from them, the principal normalised eigenvectors were derived to represent the orientation of fibrils based on the assumption that the diffusion of water in the direction of fibres is faster than that in the perpendicular direction (Le Bihan et al., 2001). Principal eigenvectors from a portion of the cartilage tissue corresponding to 74 columns and 17 rows (approximately 5.18 mm length, 1.19 mm depth) were then isolated to implement in the FE model.

Models
The portion of the bovine cartilage specimen with dimensions of 5.18 mm length (x direction in Fig. 1) Â 2 mm width ( y direction) Â 1.19 mm thickness (the z direction) was considered. Five biphasic material models, each of which included a different level of complexity of fibril reinforcement, were individually applied to the cartilage specimen. In Model 1, the primary orientation of the collagen fibril was derived from DT-MRI. The FE cartilage model was meshed into 74 and 17 elements in the length and thickness (depth) directions respectively and the principal eigenvectors derived at by DT-MRI (Section 2.1) were mapped to the finite elements of the model to represent the realistic primary orientation of collagen fibres. The statistical distribution of collagen fibres around the primary direction (Federico and Gasser, 2010;Federico and Herzog, 2008c;Wu et al., 2016) was not considered in this study. In the width direction, 32 elements were used, and it was assumed the orientation of collagen fibres remained constant through the width. In model 2, the inhomogeneous distribution of the collagen fibrils was represented by the widely-used three-layered model. The primary direction of the fibrils was assumed to be in the x direction in the SZ (parallel to the articular surface) and in the z direction in the DZ (perpendicular to the articular surface). In the SZ, the two orthogonal directions were in the y and z direction while in the DZ, they were in the x and y directions. The fibrils in the MZ were equally distributed along the three perpendicular axes, which were aligned with the global x, y and z directions. The thicknesses of the SZ, MZ and DZ were estimated from the DT-MRI images. Models 3-5 were control models with different levels of fibril-reinforcement. Due to the complexity of considering the inhomogeneous distribution of collagen fibrils, uniformly distributed collagen fibrils are also widely used in FE models to study cartilage mechanics. Therefore, in Model 3, the primary direction of collagen fibrils was assumed to be aligned with the x direction throughout the model. In Model 4, only the SZ was fibril-reinforced, with the primary fibril orientation in the x direction, and the remaining zones were isotropic. In model 5, fibril reinforcement was not included. These models are summarised in Table 1.
The bottom of the cartilage was fully fixed to simulate an ideal bond to the bone. The peripheral surfaces of the cartilage were free-draining. A rigid impermeable plate (3.9 mm radius) was assumed to be in frictionless contact with the articular surface of the cartilage and a ramp load (0.6 N) was applied to the plate ( Fig. 1) in 2 s and then kept constant for 1200 s. It should be noted that physiological loading and unloading conditions of the articular cartilage are much faster (Kutzner et al., 2010). However, the creep loading condition applied here is helpful for identifying the mechanical behaviour of cartilage (Lu and Mow, 2008).

Materials
Since the main focus of this study was on the effect of the inhomogeneous distribution of collagen fibrils, to simplify the models, the dependence of the permeability on strain (Holmes and Mow, 1990) and collagen direction (Ateshian and Weiss, 2010;Federico and Herzog, 2008a) and the inhomogeneity of the compressive modulus (Schinagl et al., 1997) were not considered. The solid phase of the cartilage was modelled as a mixture of elastic non-fibrillar matrix and collagen fibres such that the calculated stress tensor was the sum of the stress tensors for the two constituents. The equilibrium compressive modulus, Poison's ratio and permeability were taken as 0.47 MPa, 0.24 and 0.0014 mm 4 / Ns (Athanasiou et al., 1991), respectively. The effective tensile moduli of the fibrillar component were taken as 6 MPa in the primary direction, which is within the range of tensile modulus of the superficial layer tested in the literature (Verteramo and Seedhom, 2004), and one third of this value in the other two directions orthogonal to the primary direction (Kazemi et al., 2011;Li et al., 2009). In the three-layered model (Model 2), the tensile moduli of collagen fibrils in the primary direction at the deep zone was assumed to be the same as that of the superficial layer, while in the equally distributed layer, it was taken as half of that in the superficial layer (Dabiri and Li, 2013;Li et al., 2009) in all three directions. The material properties of the five models are summarised in Table 1.
The five models were all solved using both Abaqus (Version 6.13, Dassault Systems Simulia Corp, USA) and FEBio (Version 1.6.0, www.febio.org) (Maas et al., 2012;Meng et al., 2013). In Abaqus, an isotropic linear elastic material was used for the non-fibrillar matrix of the cartilage.
The formulation developed by Li et al. (2009) was implemented to model the mechanical behaviour of the collagen fibrils (Eq. (1)).  Table 1 Material properties adopted for the five models investigated in this study.

Primary fibril orientation
Effective tensile modulus of the fibrillar components (Verteramo andSeedhom, 2004, Li et al., 2009) Compressive properties and permeability (Athanasiou et al., 1991) Model where E X f is the Young's modulus of the fibres in the primary direction, ε X is the strain of the fibres, E X 0 and ε E X are elastic constants, and X represents the primary direction of the fibre at the given point. The same function was also used for the orthogonal Y and Z directions. To provide an approximately equivalent tensile modulus for the fibrils to those of the FEBio model, the dependence of tensile moduli on strain was not considered in this study, i.e.
In FEBio, the neo-Hookean material model was used for the non-fibrillar matrix. The strain energy function of the neo-Hookean material is defined as (Maas et al., 2012(Maas et al., , 2013 ( ) where I 1 is the first invariants of the right Cauchy-Green deformation tensor, J is the determinant of the deformation gradient tensor, λ and μ are the Lamé parameters, related to Young's modulus E and Poisson's ratio ν as follows A previous study showed that the neo-Hookean material in FEBio and the linear elastic material in Abaqus predicted almost identical stress values for applied strains up to 20% in a uniaxial analysis (Meng et al., 2013). The strain energy function for each fibril bundle followed an exponential power law (Ateshian et al., 2009;Maas et al., 2013;Meng et al., 2014) where β is the power of exponential argument, I n is the square of the fibre stretch and ξ is the measure of the fibre tensile modulus.
For the case of β¼2 as adopted in this study, the elasticity of the fibre at the strain origin is 4ξ (Maas et al., 2013;Meng et al., 2014).

Solutions
In Abaqus, the cartilage was discretized with C3D8RP (8-node trilinear displacement and pore pressure, reduced integration) elements. The soils, consolidation analysis procedure was used since cartilage was considered as a fluid-filled porous medium. In order to reduce computational time, an appropriate value (less than 6% of the maximum pore pressure (Goldsmith et al., 1995)) was specified for UTOL, which is an automatic time incrementation switch. The node to surface contact discretization and finite sliding tracking approach were used for the models. The numerical procedure and Fortran subroutine developed by Li et al. (Li et al., 2009) were used to model the reinforcement of the fibrils.
In FEBio, the hex8 (8-node hexahedral) elements were used to discretize the modelled cartilage sample. The biphasic analysis step was used to solve the models. The sliding2 implementation was used for the contact interface. The penalty method was used to enforce the contact constraints. The auto-penalty was applied for all models to calculate a suitable initial value for the penalty factor.
Fluid load support on the articular surface, which is the ratio of the load supported by fluid pressure to the total load, was analysed in this study because it plays an important role in cartilage function in terms of reducing friction. The other mechanical outputs analysed in this study included the 1st and 3rd principal strain and stress values and the maximum shear strain and stress, since they are potential mechanical criteria for cartilage damage.

Results
The data associated with this paper are openly available from the University of Leeds Data Repository (Meng et al., 2016).
The angle of the principal diffusion eigenvector at a pixel with respect to the normal to the articular surface was referred to as the fibre angle in this study. The fibre angle through the bovine cartilage sample is plotted against depth measured from the articular surface in Fig. 2. As stated above, the specimen had 17 pixels in the depth direction and each depth corresponded to 74 pixels in the length direction. Each point in Fig. 2 is the average of all the 74 pixels at a given depth. The error bars are the standard deviation in the fibre angle of all the 74 pixels at that depth. The average angles of the principal eigenvectors of the diffusion tensors appear to be at around 75°to the normal near the articular surface while approximately 30°to the normal in the deep zone. In the middle zone, the average fibre angle continuously varied from the larger value at the superficial zone to the smaller value at the deep zone. This increasing trend from the deep zone to the articular surface is consistent with previous studies (de Visser et al., 2008b;Meder et al., 2006).
The fibre angle distribution in each of the zones is further illustrated in Fig. 3, where the fibre orientations are shown as a proportion of fibres within each zone. In the superficial zone, 85% of the fibres were orientated at 70-90°to the normal to the articular surface. The fibres at the middle zone appeared to be equally distributed within the range of 0-90°without a preferable direction. In the deep zone, 70% of the pixels had a fibre angle less than 40°. Consistent with the previous studies (Filidoro et al., 2005;Meder et al., 2006), the zone just above the tidemark also appeared to have no preferable direction.
The fluid pressure distributions predicted by Abaqus and FEBio for the five models are shown in Fig. 4. Generally, FEBio and Abaqus produced very similar results for the five models, providing confidence that the five models in this study were all mathematically correctly solved.
The interstitial fluid load support at the articular surface of the five models is shown in Fig. 5. The fluid load support of the threelayered model (Model 2) and the uniform fibril model (Model 3) was very close to that of the DT-MRI model (Model 1) during the whole creep period, all of which reduced from 80% to less than 5% after 1200 s of creep. When the load was just applied (2 s), the fluid load support of Model 4 was 12.5% lower than Models 1-3 but 75% higher than Model 5. Consistent with the fluid load support, the magnitude of the fluid pressure of Models 2 and 3 was generally close to that of Model 1, while the fluid pressure at the articular surface of Model 4 was considerably smaller than that of Models 1-3 but larger than that of Model 5 (Figs. 4 and 5).
Although the peak magnitude of the 1st principal (maximum tensile) strain of Models 1-3 occurred at the cartilage-bone interface, Model 1 produced considerably higher values at the articular surface (almost at the same level as the peak magnitude) than Models 2 and 3 (Fig. 6). Moreover, due to the different levels of fibril reinforcement, Model 4 predicted higher values for the 1st principal strain in the MZ; and the 1st principal strain of the articular surface of Model 5 was even higher than that at the cartilage-bone interface (Fig. 6). The peak value of the 1st principal stress of Model 1 was considerably higher than that of the other models (Fig. 7). Moreover, the location of the peak 1st principal stress of the five models was also different: in the DT-MRI model, the peak occurred on the articular surface, while for the other four models, it occurred on the edge between the articular surface and the peripheral surfaces (Fig. 7).
The distribution patterns of the 3rd principal (maximum compressive) strain (Fig. 8) and stress ( Fig. 9) of Model 1 were not different from those of the other four models. For all the five models, the peak magnitudes for the 3rd principal strain and stress occurred at the cartilage-bone interface. The 3rd principal strain of Model 4 on the articular surface was considerably lower than that of Model 5 while higher than those of Models 2 and 3 Fig. 3. Fibril angle distribution within each zone. Using the results presented in Fig. 2, the superficial zone was defined as the top two layers, the middle zone from the 3rd to the 8th layer, the deep zone from the 9th to the 15th layer, while the 16th and 17th layers were considered as a separate zone.       ( Fig. 8(k)). This comparison indicated that the SZ, MZ and DZ fibrils all considerably contributed to the effective compressive stiffness of cartilage.
The peak magnitudes of the maximum shear strain of the five models all occurred at the cartilage-bone interface. The peak value of the maximum shear strain of Model 3 was considerably lower than those of the other four models (Fig. 10). At the articular surface, Models 1, 4 and 5 predicted higher maximum shear strain than Models 2 and 3 (Fig. 10). Moreover, compared with Models 1-3, Model 4 produced higher values for the maximum shear strain in the MZ. Although Models 1-4 predicted peak maximum shear stress at the articular surface, the peak magnitude predicted by Model 1 was considerably higher than those of Models 2-4 (Fig. 11). Differently to the other four models, the peak maximum shear stress of Model 5 was at the cartilage-bone interface.

Discussion
The distribution of collagen fibrils in AC is highly inhomogeneous, varying with both depth and location. Although this inhomogeneity may influence the biphasic mechanics of AC, the effects have mainly been investigated using simplified threelayered models. Therefore, the aim of this study was to investigate the effect of the realistic collagen fibril orientation on the biphasic mechanics of articular cartilage.
Most previous experimental studies that have characterised the fibrillar structure through the depth of the articular layer have only reported the dominant direction of the fibre orientation and there is limited data on the fibre angular distribution in each zone. The analysis in this study on the distribution of cartilage fibre orientations in each zone therefore provides additional information on the fibre angles through the structure. The collagen orientations obtained from DT-MRI images in this study were consistent with previous studies. Using the DT-MRI techniques, Meder et al. (2006) and de de Visser et al. (2008b) reported that the principal component of the diffusion tensor near the articular surface appeared to be oriented at around 70-75°to the normal to the articular surface while in the deep zone it appeared to be at approximately 20°, agreeing with the orientations found in this study (Fig. 2). Moreover, the collagen angle profile obtained in this study was also close to the findings measured by other imaging techniques. For example, Xia et al. (2001) measured the variation in collagen angles with depth using a more established technique, polarised light microscopy, and reported a collagen angle profile that was close to the results in this study. A similar trend was also reported by Mollenhauer et al. (2003) using X-ray diffraction. These consistencies suggests that the collagen orientations used for Model 1 were reasonably realistic. There were also some differences between this work and previous studies. For example, the dependence of the standard deviations of collagen angles on depth was not as obvious as in a previous study (de Visser et al., 2008a), where there were larger variances in the middle zone than at the peripheries. The differences might be caused by different tissue species, exact anatomical location of the sample, age, level of degeneration, freeze-thawing and storage conditions. Indeed, even for the same type of subject (bovine), de Visser et al. found a considerable degree of variability in the depth profiles of fibril angles (de Visser et al., 2008a).
Besides the above consistencies of the collagen fibre orientations obtained in this study with previous studies, the validity of the results and conclusions of this study was also supported by two other findings. First, FEBio and Abaqus predicted similar results for all the five models (Fig. 4). This agreement between the results obtained by two different software packages provides confidence that the solutions of this study are computationally correct. Secondly, the conclusions drawn from the comparison between the control models are consistent with previous studies. For example, compared with Model 5, the inclusion of the fibril reinforcement in Models 2-4 constrained the lateral deformation and reduced the tensile strain at the articular surface (Fig. 6). Moreover, the comparison between Models 4 and 5 showed that the SZ fibrils not only reduced the tensile strain on the articular surface, but also increased the compressive stiffness (Fig. 8), fluid load support (Fig. 5) and tensile stress (Fig. 7) on the articular surface, consistent with the previous studies (Korhonen et al., 2002;Owen and Wayne, 2011).
The comparison of the fluid load support between Models 1 and 2 showed that the anisotropic distribution of collagen fibrils itself does not appear to contribute considerably to the interstitial fluid load support: the three-layered model produced similar fluid load support to the DT-MRI model (Fig. 5). Therefore, the threelayered model can be used to predict the fluid load support. It should be noted that, in this study, the three-layered model (Model 2) predicted similar fluid load support to the model where the fibre alignment was the same through all layers (Model 3). This result does not seem consistent with the previous studies, in which fluid load support of the three-layered model was higher than that of the uniform fibril model (Dabiri and Li, 2013;Krishnan et al., 2003). This is because the depth-dependent variation in the compressive stiffness of the non-fibrillar matrix, which contributes to the fluid load support together with the tensile stiffness of the fibrils , was not considered in this study. Although the location of the peak 1st principal stress and maximum shear stress of the DT-MRI model (Model 1), the threelayered model (Model 2), and the uniform fibril model (Model 3) all occurred at the articular surface, the peak magnitudes of the DT-MRI model were substantially larger than Models 2 and 3 ( Fig. 7(k) and Fig. 11(k)). Moreover, although the peak 1st principal strain and maximum shear strain values of Models 1-3 all occurred at the bone-cartilage interface, Model 1 produced considerably larger magnitudes for these strains than Models 2 and 3 on the articular surface (almost similar to the peak values) (Figs. 6 and 10). Such differences in the stress and strain distributions between Model 1 and Models 2-3 were caused by the disorder in the orientation of the collagen fibre bundles. The maximal shear stresses (Atkinson et al., 1998;Eberhardt et al., 1991), excessive tensile stresses (Kelly and O'Connor, 1996), excessive principal strains (Wilson et al., 2006), and shear strain (Wilson et al., 2006) have all been suggested as damage criteria. Moreover, the typical method to estimate a damage mechanism is to compare the location of damage with the peak magnitude of a mechanical parameter (Atkinson et al., 1998;Wilson et al., 2006). Therefore, the above differences in the peak magnitudes and their locations for these mechanical parameters are of particular importance because they may affect the prediction of localised mechanically induced cartilage damage.
The comparison between Models 2 and 4 showed that it is important to consider the integrity of collagen fibrils for the accurate prediction of cartilage mechanics. The importance of the superficial layer to the function of cartilage has been extensively studied (Bevill et al., 2010;Hosseini et al., 2014;Wayne, 2006, 2011). The functions of MZ and DZ fibrils have been less extensively investigated, although they may also play important roles in cartilage mechanics (Shirazi and Shirazi-Adl, 2008). This study showed that MZ and DZ fibrils also considerably contribute to the fluid pressurization by constraining the lateral deformation of cartilage (Figs. 4 and 5). Moreover, ignoring the MZ and DZ fibrils also caused the increase in the 1st principal strain and the maximum shear strain in the MZ (Figs. 6 and 10) and the decrease in the effective compressive stiffness (Fig. 8). Therefore, it is important to achieve the integrity of collagen fibrils in engineered cartilage.
There are limitations in this study. First, the DT-MRI data was acquired from a single 2 mm slice. The DT-MRI model was created based on the data of this slice, assuming that all fibre orientations in the width direction-the direction perpendicular to the slicefollowed the same orientation. Moreover, it has been reported that the orientation of collagen fibres on the planes parallel and perpendicular to the split line may be different (Jeffery et al., 1991). However, in this study an arbitrary plane was selected for scanning. Therefore, the DT-MRI model may have underestimated the spread of the collagen fibrils in the width direction of the cartilage specimen and this limitation will be addressed in future studies.
Secondly, it has been tacitly assumed that the DT-MRI eigenvectors accurately represent the collagen orientations. It is known from previous work (de Visser et al., 2008a(de Visser et al., , 2008bFilidoro et al., 2005;Meder et al., 2006) that the principal eigenvectors do approximately reflect the trend of the fibre orientations, but a proportion of the spread of the fibre orientations (which we have considered here to be an important feature) is probably a consequence of noise or inaccuracies in the DT-MRI data which will overestimate the distribution of orientations. Moreover, a higher in-plane resolution, which may also affect the MRI results (Xia, 2007), was not adopted due to the extended acquisition time. Therefore, the effect of a higher in-plane resolution is still not clear. Thirdly, since the width of the cartilage slice was only 2 mm, there is greater diffusion in the model in this direction than in intact tissue, reducing the fluid pressurization. Indeed, the highest fluid load support in this study was 80%, whereas values found in the literature  tend to be over 90%.
In addition, collagen fibrils have important effects on the permeability since the diffusion of water molecules is greater in the direction parallel to fibrils than in directions orthogonal to the fibrils. Therefore, the permeability of cartilage is not only depth dependent but also anisotropic Herzog, 2008a, 2008b;Maroudas and Bullough, 1968). The dependence of permeability on the orientation of collagen fibrils could be considered through a more complex method in which both the eigenvectors and eigenvalues are incorporated (Pierce et al., 2013). Moreover, the permeability is also strain-dependent (Holmes and Mow, 1990) and the compressive stiffness of the non-fibrillar matrix is depthdependent (Schinagl et al., 1997). Since the focus of this study was on the effect of the inhomogeneous distribution of collagen fibrils itself, the above factors were not considered.
Despite the above limitations, this study provides more insights into cartilage mechanics: the three-layered model can be used to predict the fluid load support of cartilage; however, if the localised damage of the cartilage layer is to be simulated, then it may be necessary to incorporate a realistic orientation of the fibres.

Conclusions
The effect of the realistic orientation of collagen fibrils on the biphasic mechanics of articular cartilage was investigated in this study. It was found that the widely-used three-layered inhomogeneous model and the uniformly oriented fibril-reinforced model can produce similar fluid load support to the model including the realistic fibril orientation; however, an accurate prediction of the other mechanical parameters requires the inclusion of the realistic orientation of collagen fibrils.  d   1  1  1  1  3  3  2  2   1  1  1  1  3  3  2  2   1  1  1  1  2  2  3  3 where = a 1/ 3 1 , a 2 ¼0.9342, a 2 ¼0.3568, and where we have added a null-direction to indicate the absence of diffusion gradients (the " = b 0" images). For each non-null direction, positive and negative amplitudes were used, giving twenty gradient directions. For each of these twenty directions, four diffusion gradient values of 24.3, 34.4, 42.1 and 48.6 mT/m ( Δ = 40 ms, δ = 10 ms) were used corresponding to b-values of approximately 155, 310, 465 and 620 s/mm 2 . Three = b 0 images were also acquired, giving a total of 83 images.
However, the presence of spoiler gradients complicates and significantly affects the encoding matrix and b-values. where = a 1/ 2 4 . The spoiler gradient amplitude was fixed at 117.7 mT/m, and had a b-value of 114 s/mm 2 (Δ = 6.28 sp ms, δ = 5 sp ms). The spoiler gradients therefore modified the overall range of b-values. The minimum b-value, corresponding to the " = b 0" images was 114 s/mm 2 (i.e., only spoiler gradients), and the maximum b-value was 970 s/mm 2 .

Appendix B. Calculation of the diffusion tensor
The essentials of the DTI sequence timings are shown in Fig. A-1. Shown are the diffusion gradients (which may be applied in any direction), the spoiler gradients, and unknown constant background gradients which are assumed to be spatially linear on the scale of a pixel. All other imaging gradients occur in the time periods t 1 and t 2 except the slice-selection gradient for the re-focusing 180°pulse (not shown). This slice-selection gradient is of short duration and low amplitude and makes a negligible contribution to the calculation of the b-tensor. The other imaging gradients, being outside the diffusion-encoding period, also do not contribute to the b-tensor and can also be ignored. Thus, in order to calculate the diffusion tensor, only the effects of the three types of gradients are considered: diffusion gradients, spoilers, and a constant background gradient (as shown in Fig. A-1).
In the second-order diffusion tensor model (Basser et al., 1994), the NMR signal attenuation is given by (with Einstein summation convention on repeated indices) by negating all amplitudes either prior to the refocusing pulse or in the period afterwards. D ij is the diffusion tensor which is assumed to be time-independent. Note that, for convenience, the gyromagnetic ratio γ has been absorbed into the definition of ( )  Schematic timings of the DTI sequence (not to scale). The waveforms of the diffusion gradients and spoiler gradients are shown in relation to the RF pulses of the spin-echo sequence, together with an added background gradient (assumed linear). All other imaging gradients are not shown and occur in either the interval t 1 or t 2 with the exception of the slice-selection gradient (also not shown) of the 180°r efocusing pulse.