Non‐affine fiber kinematics in arterial mechanics: a continuum micromechanical investigation

There is growing experimental evidence for non‐affine deformations occurring in different types of fibrous soft tissues; meaning that the fiber orientations do not follow the macroscopic deformation gradient. Suitable mathematical modeling of this phenomenon is an open challenge, which we here tackle in the framework of continuum micromechanics. From a rate‐based analogon of Eshelby's inhomogeneity problem, we derive strain and spin concentration tensors relating macroscopic strain rate tensors applied to the boundaries of a Representative Volume Element (RVE), to strain rates and spins within the tissue microstructure, in particular those associated with fiber rotations due to external mechanical loading. After presenting suitable algorithms for integrating the resulting rate‐type governing equations, a first relevance check of the novel modeling approach is undertaken, by comparison of model results to recent experiments performed on the adventitia layer of rabbit carotid tissue.


INTRODUCTION
Thermodynamically sound continuum mechanics models for soft tissues, as reviewed and developed over the last decades by Holzapfel and co-workers, [20,30,40] have resulted in a thriving research community with increasing impact in biomedical engineering and biomedicine, and particularly so in cardiology. [32,33,35,47] In this context, the importance of tissue anisotropy has been more and more understood and considered, thereby aiming at a more and more resolved representation of the collagen and elastin fiber morphologies evolving during tissue deformation. Most of the corresponding contributions tackle the problem of fiber re-orientation through the use of affine transformations, i.e. the fiber movements are directly linked to the "macroscopic" strain tensor characteristics of the representative volume element associated with a piece of tissue. While this concept frequently provided very satisfactory results, in particular so in the context of mitral valve leaflet modeling, [39] we also note several experimental observations where the fibers do not follow such a deformation pattern. These observations concern the adventitia layer of carotid arteries, [36,37] but also tissues beyond the cardiac realm, such as the human liver capsule and murine skin. [34] Typically, the aforementioned deformation patterns are associated with large shear strains in the soft matrix being situated in-between the fibers; and such discrepancies between macroscopic and microscopic strains, being incompatible with affine deformation characteristics, have also been reported for tendon fascicles. [23,24,55] These interesting observations have motivated the theoretical developments described in the present paper: Here, the continuum micromechanics of fiber-matrix composites and polycrystals, which turned out as a versatile tool for predicting the material behavior of biomedical materials such as bone and bone biomaterials, [16][17][18]45,59] is consistently extended towards the representation of materials undergoing large deformations, with corresponding load-induced changes in the fiber morphology, i.e. large fiber rotations.
Accordingly, the paper is organized as follows: In Section 2, we extend the classical notions of continuum micromechanics, from (linearized) strains to strain rates, providing the basis for a continuum micromechanics theory applicable to microheterogeneous materials undergoing large deformations and large (micro-)configurational changes. Thereafter, in Section 3, we specify these developments for (rotating) fiber-matrix morphologies, employing an extended version of Eshelby's matrixinclusion problem, which provides access to concentration (or down-scaling) tensors, not only for strain rates, but also for strain rate-spin interactions. Section 4 is devoted to the integration of the aforementioned rate-type equations, by means of suitable algorithms. Section 5 provides corresponding computational results, allowing for the comparison of fiber rotations predicted from the novel micromechanical approach developed herein, with fiber rotations obeying the classically assumed affine transformations. Conclusions are finally drawn in Section 6.

Operator Name
(.) 0 quantity (.) in the initial configuration ⟨(.)⟩ average of quantity (.) over the RVE (.) transpose of tensorial quantity (.) (.) −1 inverse of tensorial quantity (.) (.)∕ partial derivative of quantity (.) with respect to time (.)∕ material derivative of quantity (.) ⋅ inner product (first-order tensor contraction) : second-order tensor contraction ⊗ dyadic product ∇ objective logarithmic derivative of det determinant operator div divergence operator macroscopic Eulerian gradient with respect to the current location vector microscopic Eulerian gradient with respect to the current location vector Variable Name second-order identity tensor ( ) fourth-order strain rate concentration tensor ∞ strain rate concentration tensor of the Eshelby problem

Representative volume element -average rules
We consider a piece of fiber-reinforced soft biological tissue (such as arterial tissue), having potentially undergone large configuration changes, and filling, in the current configuration, a so-called representative volume element (RVE). Definition of the latter rests upon the separation of scales requirement. [12,65] Mathematically, this reads as with as the characteristic length of RVE (hundreds of micrometers in the case of arterial tissue), as the size of heterogeneities within the RVE (relating to 10 micrometer fiber diameter in arterial tissue), and  as the structural length associated with the investigated mechanical problems [1] ; e.g.  = ‖ ( )‖∕‖|GRAD ( )|‖ in the case of a macroscopic stress field ( ) occurring at the organ scale of several millimeters to centimeters, with being the macroscopic location vector; and GRAD being the gradient operator at the macroscopic scale, quantifying changes arising from scanning across the organ scale. Extending the notions given by Hashin, [26] from strains to strain rates, a homogeneous ("macroscopic") strain rate tensor is prescribed, in terms of a velocity field ( ) across the boundary of the RVE. Mathematically, this reads as with labeling (microscopic) points inside the RVE and on its (current) boundary, the latter being denoted as Ω. Inside the RVE, the velocity field is continuously differentiable, inducing microscopic strain rates ( ) of the format with Ω as the (current) volume of the RVE, as the microscopic gradient, and as the transpose operator. Equations 2 and 3 directly imply a strain rate average rule analogous to the strain average rule given by Hashin, [26] reading as where ⊗ denotes the dyadic product, and as the outward unit normal vector at point of the boundary of the RVE. Furthermore, the aforementioned strain rates provoke traction forces ( ) on the boundary of the RVE. These traction forces need to be in equilibrium, regardless of the particular shape and size of the RVE, leading to whereby use of Cauchy's theorem ( ) = ( ) ⋅ ( ) has been made. The external power  exerted by the aforementioned traction forces reads as whereby use of Equation 5 was made. Equation 6 manifests that the force quantity exerting power on the macroscopic strain rates is the volume integral over the microscopic Cauchy stress, which is independent of microscopic position, and of dimension "stress times volume". This induces the existence of the macroscopic Cauchy stress in the form: Equation 7 represents the well-known stress average rule. Insertion of (7) into the principle of virtual power [22,53] yields with  as the power of internal forces. Equation 8 is an analogue to the so-called Hill's lemma, stating the equivalence of the power exerted by the macroscopic stresses on the macroscropic strain rate, and the average over the power exerted by the microscropic stresses on the microscopic strain rates. Mathematically, this reads as:

Hypoelastic material behavior at the microscopic scale
The matter within the RVE behaves hypo-elastically, [58] i.e. (micro-)strain rates and (objective) stress rates ∇ are related by location-dependent elasticity tensors ( ), reading as Among the many objective stress rates documented in the open literature, we choose the logarithmic stress rate proposed by Xiao and colleagues, [61][62][63][64] since it guarantees simultaneous objectivity with respect to both stress and power, and may also provide, if desired, a direct path towards an hyperelastic analogue, as outlined in [63]. The corresponding logarithmic rate ∇ of a quantity is defined by: with ∕ denoting the material derivative of quantity , and with ( ) as the local material rotation tensor (also called logarithmic spin tensor), which is defined by: In Equation 12, ( ) stands for the spin tensor, being defined as the antisymmetric part of the velocity gradient and ( ) is given by (see [62] for full proof of this expression): In expression (14), all quantities are related to the point within the representative volume element; is the second-order microscopic left Cauchy-Green tensor with eigenvalues 1 , 2 , and 3 . is defined as a function of the microscopic transformation gradient , namely Furthermore, the coefficient in Equation 14 is defined through and the coefficients are defined through The question arises of how to relate the microscopic strain rate and spin fields to the macroscopically imposed loading. The equations governing the mechanical problem of the RVE, being related to equilibrium (7), compatibility (4), boundary conditions (2), and material behavior (10), are all linear with respect to the prescribed strain rate. Consequently, the microscopic strain rate and the microscopic spin tensor are linearly related to the prescribed boundary conditions: with ( ) and ( ) as the microscopic fields of strain rate and of spin tensor, respectively, ( ) as the fourth-order strain rate concentration tensor, and ( ) as the fourth-order strain rate-to-spin concentration tensor. In order to arrive at explicit expressions for  and , we specify the morphology inside the RVE as one of rotating fibers being embedded into a soft, compliant matrix; and make extended use of Eshelby's matrix inclusion problem, as described next.

ESHELBY PROBLEM-BASED MICROMECHANICS OF LARGE FIBER ROTATIONS IN SOFT MATRICES
In this section, we derive analytical expressions for the strain rate and spin concentration tensors ( ) and ( ); for an RVE hosting several fiber phases which are embedded into a soft matrix phase, see Figure 1(a). Starting point for these derivations is the matrix-inclusion problem of Eshelby. [13] Based on an extension of Eshelby's solutions towards the evaluation of spin tensors, we propose a new variant of the so-called Mori-Tanaka scheme. [2,44] F I G U R E 1 (a) Representative Volume Element of soft tissue subjected to homogeneous strain rate ; (b) definition of the fiber orientation through Euler angles and

RVE of rotating fiber-reinforced soft tissue
Adopting the classical strategy of continuum micromechanics, [65] we consider the case where the function ( ) in Equation 10 can be split into several uniform stiffness fields, called material phases. More precisely, we introduce one matrix phase and fiber phases. The latter are characterized by volume fractions , , = 1, … , , and stiffness tensors ; while the matrix phase exhibits stiffness and volume fraction = 1 − ∑ , . The orientation of the -th fiber phase is defined through Euler angles and , as seen in Figure 1(b). This problem definition naturally induces the existence of phase-specific strain rate and strain rate-to-spin concentration tensors, i.e. phase specific-tensors , , ℝ , , as well as need to be determined, rather than tensor fields ( ) and ( ). The phase-specific concentration tensors relate the macroscopic strain rate tensor to the strain rate and spin tensors averaged over the phases, so that Mathematical expressions for , , ℝ , , and can be derived in a particularly elegant way, through identification of the fiber phases with prolate spheroidal inclusions (or inhomogeneities) occurring in the so-called Eshelby problem, which is described next.

Eshelby's inclusion problem -definition of the rotation operator
We here re-develop the original Eshelby problem in terms of strain rates rather than strains, in order to provide the basis for a large strain micromechanics theory. The latter also encompasses the consideration of configurational changes, and therefore we are not only interested in the strain state of the inclusion, but also in its rotation or spin. Accordingly, we consider an infinitely extended, homogeneous matrix with elastic stiffness ℂ 0 , a subvolume of which, being denoted by , undergoes a uniform eigenstrain rate , while the rest of the matrix is free from any type of loading, see Figure 2(a).
This situation provokes a uniform velocity gradient field in the inclusion, reading as with ℎ as the strain rate-to-velocity gradient tensor of the Eshelby problem, depending on the shape and the orientation of the ellipsoidal inclusion, and on the stiffness tensor ℂ 0 of the matrix. For infinitely long, prolate spheroids, representing rotating fibers in the current context, and embedded in an isotropic matrix, ℎ takes the following format in terms of non-zero components with respect to the ( , , ) base frame of Figure 1 with 0 as the hypoelastic Poisson's ratio of the matrix, see the Appendix for more details on the mathematical derivations of (21). As a direct consequence of the uniformity of the velocity gradient within the inclusion, both the strain rate tensors ( ) and the spin tensors ( ) are uniform within the inclusion and linearly related to the prescribed eigenstrain rate , through: with ℎ and ℝ ℎ as the fourth-order Eshelby tensors, the first one being related to inclusion strain rates and the second one to inclusion rotations Accordingly, the tensors ℎ and ℝ ℎ depend on the shape and orientation of the ellipsoidal inclusion, as well as the stiffness tensor of the matrix ℂ 0 . Outside the inclusion, the strain rates and spins are not uniform, and in this context, particularly strong variations in strain rates and spins are encountered close to the inclusion. This inclusion problem can then be extended towards the inhomogeneity problem, related to the case of an infinite homogeneous material with stiffness ℂ 0 , in which a small ellipsoidal inhomogeneity with stiffness ℂ is embedded, while the infinite matrix is subjected to a uniform strain rate 0 at infinity, see Figure 2(b). In fact, total identity of the stress and strain rate fields prevailing in this new problem, with those of the original inclusion problem, can be gained through setting of As a consequence, the associated strain rate and spin fields remain uniform within the ellipsoidal inclusion, and their values are linearly related to the applied loading: whereby ℝ ℎ is still defined through Equation 24, and the Hill tensor ℙ reads as The fourth-order tensor ℝ ∞ is the operator which linearly relates the microscopic spin tensor to the applied macroscopic strain rate 0 , and it is given by: whilst the fourth-order tensor ∞ is the operator linearly relating the microscopic strain rate tensor to the applied macroscopic strain rate 0 : The strain rate 0 , imposed at infinity, to the Eshelby problem of Figure 2(b), needs to be related to the strain rate subjected to the boundary of the RVE of the fiber-reinforced soft tissue. This is tackled next.

Mori-Tanaka scheme for matrix with embedded, rotating fibers
In order to relate the auxiliary quantity 0 acting on the Eshelby problem, to the strain rate subjected to the RVE of Figure (1), we adopt the strategy outlined in, [2,44,65] and associate with each fiber phase, an Eshelby problem with an inclusion exhibiting the fiber stiffness tensor and shape, such that the uniform inclusion strain rates and spins are identified with those of the corresponding fiber phase, and an auxiliary matrix exhibiting the matrix stiffness tensor. This implies simplification of the strain rate average rule (4) to and insertion of (26) 1 specified for all phases, into (30) yields a relation between 0 and , reading as Insertion of (31) into (26) yields the sought-after strain rate and strain rate-to-spin concentration tensors introduced in Equation 19 as whereby ∞ , and ℝ ℎ , follow from Equations 29 and 28, when identifying the (geometrical and mechanical) properties of fiber phase , as those of the inclusion ; and ∞ follows from Equation 29, when identifying the (geometrical and mechanical) properties of matrix phase, as those of the inclusion , i.e. ∞ = , and does not occur anymore in Equation 32.

TEMPORAL INTEGRATION OF GOVERNING EQUATIONS -INCREMENTAL ALGORITHM
The homogenization scheme of Sections 2 and 3 is defined in rate form; hence, temporal integration is required for predicting the material behavior over a longer time period, such as that needed for quasi-statically loading an RVE of the considered material. We here consider displacement-driven loading scenarios, i.e. we prescribe a macroscopic deformation gradient history = ( ), referring to stress and strain-free initial configuration, i.e. ( = 0) = , with as the second-order identity tensor. The relation between this deformation gradient and the Eulerian strain rate tensor can be found in any pertinent textbook on continuum mechanics [29,53] with the temporal derivative of the Green-Lagrange strain tensor reading as is evaluated at a series of discrete time instants , with = 1, … , , which are all separated by a time interval The corresponding values of are denoted as = ( ). Then, starting from = 1 , the following sequence of algorithmic steps is performed: The Eulerian strain rates are concentrated into the fiber phases, in terms of both phase-specific strain rates and spins, ⋅ , and assuming, in the sense of a forward integration scheme, the spin tensor , and the strain rate , to be constant between time instants and +1 , the base vector associated to fiber phase at time instant +1 follows as whereby +1 , is normalized. These updated base vectors give access to updated Euler angles +1 and +1 , and hence, to updated concentration tensors +1 , and ℝ +1 , . Update of the phase-specific stresses requires knowledge of the logarithmic rate occurring in Equation 10, and hence of the logarithmic spin rate occurring in Equations 11, 12, and 14. The latter follows from the phase-specific (microscopic) deformation gradient , starting at identity, ∀ ∈ {1, … , , }, ( = 1 ) = , and evolving according to which implies the following updating rule in the context of a forward integration scheme, Combination of (39) with (35), (13), and (3) yields with = + ℝ . With the help of (15), the deformation gradient +1 gives access to the left Cauchy-Green tensor +1 , and the left Cauchy-Green tensor at time , via (12), (14), (16), and (17), yields , , . This completes the collection of updated microscopic variables needed for updating of the microscopic stresses, according to the evaluation of Equations 10 and 11 at time instant , yielding Integration of the material derivative of the Cauchy stress ( ∕ ) will be based on its relation with the second Piola-Kirchhoff stress , see e.g. the book of Salencon [53] with = det as the Jacobian determinant. As is defined with respect to the non-moving reference frame, the partial and material derivatives of the second Piola-Kirchhoff tensor are identical, = and combination of (43) with (42) yields The partial derivative of with respect to time is then approximated in the context of a forward integration scheme, yielding, together with (42) The updated Piola-Kirchhoff stress tensor +1 is then converted, by means of (42), to the updated Cauchy stress, yielding Moreover, taking the determinant of (39) yields which allows for transformating Equation 47 into As a forward integration scheme is only suitable for small time increments Δ , it is natural to consider a first-order development of (49) with respect to Δ , which delivers the sought integrated format of the Cauchy stress as whereby it was considered that Combination of (50) with (41) yields Finally, the corresponding macroscopic stress tensor follows from the stress average rule, in the format: This concludes the stress update following from prescription of a constant strain rate over a time interval Δ .

APPLICATION TO FIBER ROTATION IN ARTERIES
We now apply the aforementioned algorithm to loading scenarios tested on the adventitia layer of rabbit carotid tissue, as reported by Krasny et al. [36] More precisely, we consider tension-inflation experiments, related to the following format of the macroscopic deformation gradient, with as a constant axial deformation rate, amounting to 0.2∕min, and = 0.65 as the average of the ratio between circumferentially and axially measured stretches and . [36] The latter are the eigenvalues of the right stretch tensor = √ ⋅ , and they read as In a very first approximation, the collagen fiber distribution, as quantified from multiphoton microsocopy, [36] is represented by four fiber families with a volume fraction of ,1 = ,2 = ,3 = ,4 = 4 = 1− 4 = 8.75%. The four fibers exhibit different orientations, defined through the following pairs of Euler angles: ( 1 , 1 ); ( 1 , 2 ); ( 2 , 1 ); ( 2 , 2 ). Thereby, 1 = 90 • and 2 = 270 • , since the collagen fibers lie in the plane spanned by the circumferential and the axial directions. We consider three different arterial samples (labelled by superscript ( ) ) with the following co-latitudinal angles ( ) ,0 : (1) 2,0 = 59.5 • 2,0 = 81.5 •  [36] (bars indicate standard deviation around mean value) modulus of 500 MPa, which appears as tangent modulus in the uniaxial stress-strain experiments of Sasaki and Odajima, [54] and to a hypoelastic Poisson's ratio of 0.34, a value which we adopt from conventional elastic descriptions [10,59] ), as well as by a matrix stiffness tensor of = 2 + 3 , with hypoelastic bulk modulus = 0.83 kPa and hypoelastic shear modulus = 0.38 kPa (corresponding to a very low hypoelastic modulus of 1 kPa, which is the tangent modulus of cells reported by Stricker et al. [57] ).
The results are documented in Figure 3: It becomes obvious that, despite the morphological simplifications assigned to the RVE (when compared to the "real" microstructure seen in the microscope), the rotation characteristic is satisfactorily represented by our new micromechanical model. This may be compared to results from the affine transformation assumption, i.e. to angle predictions of the format [3,38] = tan −1 It turns out that the affine transformation underestimates the fiber rotations observed in the adventitia layer of the rabbit carotid tissue. Accordingly, the shear deformations in the matrix phase need to be larger than the macroscopically prescribed ones -a model feature which is not compatible with the affine assumption.

DISCUSSION AND CONCLUSIONS
To the best knowledge of the authors, this is the first theoretical approach for non-affine fiber transformations in the setting of Eshelby-problem based continuum micromechanics. Naturally, it considers large morphological changes associated with large strains subjected to the investigated RVEs, and in this context, it constitutes an interesting hypo-elastic complement to the collection of large strain hyper-elastic formulations developed in the framework of homogenization theory. [41,42] In this context, the originality of our approach relates to the use of hypoelasticity in a micromechanical framework, while we are obviously not the first ones to advocate the use of hypoelasticity per se, as a suitable tool for soft tissue biomechanics. In this context, we may refer to the interesting contributions of Freed [14,15] who motivates the use of hypoelastic formulations on historical as well as experimental grounds. His historical argument refers to Y. C. Fung, the "father of biomechanics", who in 1967 described the nonlinear behavior of a rabbit mesentary by a first-order differential equation, [19] which can be readily recast, via the chain rule, into a hypoelastic equation. [14] Freed's experimental argument for hypo-elasticity builds on the measurements of Criscione et al., [8,9] which show the path-dependency of stress-strain relations in soft biomembranes. Such a path dependence cannot be accounted for by a hyperelastic approach, where the stress depends solely on the current total strain, irrespective of how the latter was accumulated. Hence, hypoelasticity appears as a sometimes necessary generalization of hyperelasticity, as is known from the famous 1955 contributions of Truesdell and Noll. [46,58] This generalization, however, has to be done with care, so as to avoid the violation of fundamental mechanical principles, such as objectivity and thermodynamic consistency.

Comparison of different objective stress rates
Objectivity refers to the requirement that temporal derivatives of stress eigenvalues must remain invariant under a rigid body rotation. Hence, the temporal derivatives of arbitrary Cauchy stress components need to be complemented by suitable rotation elements; and several suggestions on such complements have been proposed in the open literature, the most popular ones being the Truesdell rate, the Jaumann rate, and the Oldroyd rate, Rather than the aforementioned choices, we prefer the logarithmic rate of Xiao and colleagues, [61] as it is the only measure which provides, at the same time, objectivity of both stress and power rates. Still, the question might arise if our specific choice, based on a strong theoretical argument, would have a remarkable effect on actual simulation results. Repetition of the computations described in Section 5, while replacing the logarithmic rates (11)- (14) by each of the rates (58)-(60), shows that this is actually not the case, see Figure 4. At the same time, Figure 4(a) shows a non-linear stress-strain behavior which is classically expected for soft tissues. However, the model predictions appear as too stiff when compared to experimental observations, and we think that this open challenge might probably be overcome through a multiscale approach where the currently introduced fibers are made up of yet smaller fibers. However, this is clearly beyond the scope of the present manuscript.

Thermodynamics basis for hypoelasticity
While the term "elasticity" has undergone, throughout history, a multitude of meanings, even in the field of mechanics, [49] the most appropriate and general (mechanical) meaning of "elasticity" refers to the non-dissipative behavior of bodies under deformation. The latter is normally checked through the Clausius-Duhem expression for the dissipation ; reading for isothermal states as [50][51][52] F I G U R E 4 Stress-stretch response (left) and fiber inclination evolution (right) for different choices of the objective stress rate with the current mass density and with the Helmholtz free energy (per unit mass). The Helmholtz free energy is a state function, standardly depending on a strain-like variable and on temperature, and possibly also on other internal variables. When using the Green-Lagrange strain tensor as state variable, Equation 61 is conveniently expressed in terms of the reference configuration, reading as, see e.g., the book of Salencon [53]  = ∶̇− 0 ( ( )) = 0 with 0 as the initial mass density and with as the second Piola-Kirchhoff stress tensor; yielding the classical hyperelastic state equation Under these conditions, can also be called a "strain energy density". However, requirement (61) is met by a class of material behaviors which are much richer and more diverse than those cast in the hyperelastic formalism of strain energy density, as it was impressively shown by Rajagopal and Srinivasa. [50][51][52] In this context, it is particularly interesting to use the Gibbs potential  , which is related to the Helmholtz potential via [5,50] Insertion of (64) into the dissipation Equation 61 yields implying the following expression for the strain rate As ( ∕ ) is not objective, i.e. still depends on the motions of the observer, Equation 65 needs to be expanded by power-free or "gyroscopic" terms ∶ 2  ∶ (− ⋅ + ⋅ ) = 0 with as a spin tensor fulfilling −1 = and tr = 0. The extended format of Equation 65 then reads as yielding a thermodynamically based, objective hypo-elastic law in the format: Specification of (68) for = and for  = 1 2 ∶ −1 ∶ readily delivers the hypoelastic phase behavior of Equation 10. Hence, our choice for a hypoelastic phase behavior fulfills all the requirements concerning objectivity and thermodynamics. As compared to the hyperelastic setting, the free energy associated with hypoelasticity is "freed" from the reference configuration, and therefore "lives fully in the here and now".

Model philosophy
In accordance with the philosophical principle of Occam's razor ("the simplest explanation is best", as stated by Cover and Thomas [7] ), or Popper's quest for identification of the per se least probable, and hence, most informative theory capturing a set of observations, [48] our intention is to propose a model built on the very minimum of microstructural information, which is necessary to deliver reliable predictions of macroscopic material behavior. In addition, this microstructural information needs to be accessible by experiments; and in the latter context, the orientation of fibers appears to be far more easily accessible than the network connectivity. This was the primary motivation for not introducing the latter into our present version of the model -and we also note that such an introduction is not straightforward within the currently employed model framework. From a pragmatic as well as historical viewpoint, we note that consideration of fiber orientation (without explicit account of the network connectivity) prolongs a well-established tradition in soft tissue modeling, as pursued by Holzapfel et al., [21,31,32] among others.

Anisotropy and Eshelby's matrix-inclusion problem
Eshelby's problem relates to the behavior of an elastic inhomogeneity (which does not have to behave necessarily isotropically) when embedded in an isotropically elastic, infinitely extended matrix, which is subjected to (linearized) homogeneous strains at infinity. This (actually fictitious) infinite matrix needs to be clearly distinguished from the homogenized material within the RVE: even in case the former is isotropic, the latter may well be anisotropic. Actually, in the present case, both the fiber phases and the matrix phase behave isotropically [therefore, we present, at the phase level, values for two hypoelastic constants associated with hypoelastic isotropy, namely for the hypoelastic modulus and the hypoelastic Poisson's ratio, see the text below Equation 56], while the homogenized RVE behaves anisotropically. This anisotropy stems solely from the shape and the orientation of the fiber phases within the RVE. The latter situation is often encountered with biological tissues, not only in the example of the present paper, but e.g. also in bone. [28] As regards phase strains, it was made clear by Benveniste [2] and repeatedly discussed thereafter [56] that the Mori-Tanaka scheme identifies the average strains in the matrix phase, as those prescribed at infinity of the aforementioned fictitious matrix of the Eshelby problem.
Due to its extremely high computational efficiency, Benveniste's method of 1987 [2] has remained, over the last thirty years, the perhaps most successful and popular approach for concentration tensor estimation (as might be inferred by extremely high citation rates). However, many alternative methods qualify for concentration tensor estimation (as they would for our newly introduced rotation tensor ℝ ). Potential choice of such methods would also depend on the problems to be investigated. While a Mori-Tanaka scheme is a natural choice for matrix-inclusion morphologies, polymer-chain microstructures might be more suitably represented by Miehe's statistical multi-chain model. [43]

Affine fiber kinematics and experimental data
When referring to the limits of affine transformation modeling in the Introduction, we referred, among others, to reports on experimental studies, where the authors did not explicitly mention the word "affinity" -and this includes the papers of Screen et al. [55] In this context, it may be appropriate to review the meaning the term "affinity" in the context of kinematics and deformation of microstructures. Some authors, such as Wen et al. [60] when testing polymer hydrogels, use the term in a very broad sense, stating that an affine deformation relates to "local strain in a sample" which "is identical everywhere and equal to the macroscopic strain". This very general definition would qualify any soft tissue as non-affinely deforming, because at some level within the microstructure or nanostructure, some deviations from perfect strain homogeneity would naturally occur anyway. Thus, the question arises at which particular length scale this would actually be the case. Accordingly, the aforementioned authors do speak about "a length scale above which the gels deform affinely". Our present contribution is more focused, and we concentrate exclusively on affine versus non-affine fiber kinematics. Different experimental techniques have been employed in order to quantify the fiber strains and rotations, and those of adjacent material phases at different levels; and whenever these aforementioned "micro"-strains would be identical to the macrostrains (up to tensorial transformations which are solely related to a potential change of the employed base frame), then the fiber-matrix microstructure would obey "affine transformations". Screen et al., [55] by using a confocal scanning laser microscope, identify fiber and matrix strains at the level of tenocytes and tendon fibers within a fascicle, which differ both from the macrocospic (fascicle-related) strains, and from each other. Therefore, they identify fiber-related microstructural deformations inside a fascicle, which are clearly non-affine.

Comment on fiber interaction
We also regard the present approach as an interesting contribution to enlarging the canon of micromechanical models for fiberreinforced biological materials with strongly changing morphologies, often rooted in the frameworks of statistical mechanics and/or finite element methods, [11] towards semi-analytical homogenization schemes.
In this context, the original features of our approach are mainly concerned with modeling fiber-matrix interaction, rather than with a connected network of fibers, and with the extension of the classical analytical approaches based on the famous Eshelby problem, rather than with the ab initio development of rather "computational tools". In this context, our mode of modeling fiber interactions may deserve additional discussion: From Figures 1 and 2 it may seem that our model considers fiber interaction "only through the matrix phase". In particular, Figure 2 may insinuate that 0 is entirely prescribed onto the auxiliary matrix with stiffness ℂ 0 , which, in the case of the Mori-Tanaka scheme, is identified with the stiffness of the matrix phase, ℂ 0 = . In other words, a vanishing stiffness ℂ 0 → 0 could be expected to leave the fiber inclusion free of strain rates ( = = 0), even for non-zero values of 0 . However, deeper scrutiny into the structure of tensor ∞ according to Equation 26 reveals the following: In a base frame coinciding with the fiber direction (labelled "3" herein), the majority of the tensor components are zero: Thus, these components do not reflect any matrix properties. The same is true for 3333 , as 3333 = 1 holds irrespective of the matrix properties; while the remaining tensor components, i.e. 1212 , 1221 , 2112 , 2121 , 2332 , 3223 , 2323 , 3232 , 1331 , 3113 , 1313 , 3131 , depend on the components of ℂ 0 . Hence, even if ℂ 0 → 0, the normal strain rates 0,33 , which are subjected to the fictitious matrix in fiber direction, are fully transferred into the fiber phase, ,33 = 0,33 , even for a per se "non-acting", i.e. stiffness-free, matrix. Hence, 0,33 = ,33 is left as the only "interaction term" linking the macroscopic strain rate subjected to the RVE, to the strain rate occurring in the fiber phase. According to Equation 31, this remaining interaction term depends on the shape and the stiffness of all fiber phases, as well as on the stiffness of the matrix. We may thus conclude that a Mori-Tanaka scheme considering ellipsoidal inclusions with lengths ratios going to zero or infinity, respectively, i.e. considering "infinitely long" fiber inclusions, indeed reflects fiber interactions beyond those "through the matrix phase only". Consequently, our modeling approach does not so much relate to loosely dispersed fibers in a contiguous matrix, as rather to fibers being connected through cross-linkers, as it is the actual biological situation. We note in passing that a similar modeling concept was also successful in the realm of hard tissues: A Mori-Tanaka-type scheme associated to a cylindrical pore inclusion phase in a solid matrix phase turned to be suitable for predicting the mechanical behavior of the actually mutually connected branch-type longitudinal pore channels of cortical bone [6] ; as well as their inter-penetrating analogues in trabecular bone, [25] as reported in a series of papers. [27,28,59] Similar to these earlier approaches, our novel semi-analytical modeling tool surely entails additional intellectual challenges, but eventually provides a particularly efficient simulation tool, which we therefore regard as very helpful both for the fundamental understanding of certain types of soft biological tissues, and for the further refinement of large-scale simulations for computational medicine -thereby extending recent multiscale concepts developed for orthopaedy, see e.g. the contribution of Blanchard et al., [4] to the realm of cardiology.

APPENDIX: DETERMINATION OF THE STRAIN-TO-VELOCITY GRADIENT TENSOR
According to Eshelby's 1957 paper, [13] the velocity gradient occurring in an ellipsoidal inclusion embedded in an isotropic matrix, is linearly related to the eigenstrain rate in the inclusion, via It can be shown that the integral reduces to zero for all terms that contain an odd power of , , or . Only 21 components with respect to a base frame aligned with the axes of the ellipsoid remain, and they are equal to: