Fluid-solid interaction in the rate-dependent failure of brain tissue and biomimicking gels
Michele Terzano, Andrea Spagnoli, Daniele Dini, Antonio Elia Forte
FFluid-solid interaction in the rate-dependent failure of brain tissue andbiomimicking gels
M. Terzano a,b, ∗ , A. Spagnoli a , D. Dini b , A. E. Forte c,d a Department of Engineering and Architecture, University of Parma, Parco Area delle Scienze 181/A, 43124 Parma, Italy b Department of Mechanical Engineering, Imperial College London, Exhibition Road, London SW7 2AZ, UK c DEIB, Politecnico di Milano, Via Ponzio, 34/5 - 20133 Milano, Italy d School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts, USA
Abstract
Brain tissue is a heterogeneous material, constituted by a soft matrix filled with cerebrospinal fluid.The interactions between, and the complexity of each of these components are responsible for the non-linear rate-dependent behaviour that characterizes what is one of the most complex tissue in nature.Here, we investigate the influence of the cutting rate on the fracture properties of brain, through wirecutting experiments. We also present a model for the rate-dependent behaviour of fracture propagation insoft materials, which comprises the effects of fluid interaction through a poro-hyperelastic formulation.The method is developed in the framework of finite strain continuum mechanics, implemented in acommercial finite element code, and applied to the case of an edge-crack remotely loaded by a controlleddisplacement. Experimental and numerical results both show a toughening effect with increasing rates,which is linked to the energy dissipated by the fluid-solid interactions in the process zone ahead of thecrack.
Keywords: brain tissue; hydrogels; rate-dependent fracture; poroelasticity ∗ Corresponding author. Tel.: +39-0521-905927; fax: +39-0521-905924;e-mail: [email protected].
Preprint submitted to Journal of the Mechanical Behavior of Biomedical Materials February 23, 2021 a r X i v : . [ q - b i o . Q M ] F e b . Introduction Brain tissue is arguably one of the most complex, delicate and heterogeneous tissues of the humanbody. Its structure is characterised by two main constituents: the grey matter, which contains thenerve cell bodies, and the white matter, with a large proportion of myelinated axons (Budday et al.,2020). By a mechanical point of view, neural tissues are among the softest of all internal organs(Guimar˜aes et al., 2020), receiving protection from the skull and isolation from external actions by thecerebrospinal fluid. A large proportion of this fluid is free to move by diffusion and consolidation withinthe tissue’s solid network; as a result, the brain behaves as a soft sponge: its microstructure, albeithighly inhomogeneous, presents small pores that are saturated by fluid (Forte et al., 2017). Diffusionhas a fundamental importance for the brain function, delivering vital nutrients to the neural cells andplaying an essential role in therapies based on drug delivery (Nicholson, 2001). Besides, the motion offluid within the solid network causes volumetric shrinking and triggers consolidation effects (Franceschiniet al., 2006), which can explain various phenomena, including the onset and evolution of hydrocephalusand the brain shift during surgeries (Stastna et al., 1999; Forte et al., 2018). The interaction betweeninterstitial fluid and solid matrix provides a source of energy dissipation (Mak, 1986), which results intime-dependent behaviour, frequently observed during mechanical testing (Jin et al., 2013; Forte et al.,2017). In addition, a further source of dissipation is related to viscoelasticity, caused by intracellularinteractions between cytoplasm, nucleus and the cell membrane (Budday et al., 2017a).Mechanical models of the brain tissue at the continuum scale are usually formulated in the frameworkof finite strain mechanics, addressing the nonlinear elastic, time-dependent and hysteretical behaviour(de Rooij and Kuhl, 2016). The biphasic nature of the tissue can be captured by models derived fromthe classical theory of consolidation in soil mechanics (Biot, 1941; Franceschini et al., 2006; Forte et al.,2017), eventually coupled with large deformations (Simon, 1992; Garc´ıa and Smith, 2009; Hosseini-Faridet al., 2020). An equivalent description has been developed in the context of mixture theories (Mow et al.,1980), leading to the formulation of a consistent framework for soft porous media (Ehlers and Wagner,2015; Comellas et al., 2020). Time-dependent behaviour due to viscous effects has been described bygeneralised Maxwell models (Forte et al., 2017; Qian et al., 2018) or more refined descriptions elaboratedin the finite strain theory (Budday et al., 2017b; Haldar and Pal, 2018). However, with respect to tissuefailure, our understanding is considerably more limited. It is known that the brain tissue, as most internalorgans, does not carry significant mechanical loads; nevertheless, traumatic injuries expose the tissue todamage and fracture (El Sayed et al., 2008). Furthermore, the tissue can be perforated with catheters,needles and probes, during minimally invasive surgeries and regenerative therapeutics (Casanova et al.,2014; Ashammakhi et al., 2019; Terzano et al., 2020). Due to its high heterogeneity, failure propertiesin the brain tissue are region dependent. In the white matter, which is characterized by fibrous axonalstructures, failure occurs by tearing of fibres when the tissue is loaded above a certain threshold (Buddayet al., 2020). At the microscale, axonal injury involves a viscoelastic mechanism with stretching andsliding of microtubules, depending on the entity of the deformation (Cloots et al., 2011; Ahmadzadehet al., 2014). While the contribution of fluid-solid interaction in terms of the tissue mechanical behaviouris widely recognised, the role of fluid diffusion during failure has not been investigated. Furthermore,the flow of fluid in the brain tissue is not homogeneous (Zhan et al., 2019; Jamal et al., 2020). Dueto the different microstructure, white matter is far more permeable than grey matter, which insteadpresents densely connected networks that can entrap the fluid phase (Budday et al., 2020). Such adifference might explain the enhanced rate-dependence of white matter during compression and tensiletests, because of a faster fluid drainage (Budday et al., 2020). However, the characteristic time of fluiddraining depends on the size of the perturbed region, which also makes this contribution dependenton testing conditions (Wang and Hong, 2012a). Therefore, there is a need for investigating how rate-dependency and fluid-structure interactions affect fracture propagation in brain tissue.When a porous network is swollen by fluid, mechanical and hydraulic responses are coupled: forcesand deformations change the pressure of interstitial fluid, while pressure gradients drive fluid flow,resulting in mechanical deformation (Arroyo and Trepat, 2017). During fracture, the flow of fluid insidethe crack-tip zone might affect the surface energy required for crack initiation and propagation. Despite3carce information in the context of biological tissues, illuminating evidences come from experimentalwork on failure of hydrogels. As an example, studies on reversible gels suggest that the fracture energycan be increased by the drainage of fluid in the crack-tip zone (Baumberger et al., 2006; Naassaouiet al., 2018). In addition, polymeric gels swollen with fluid and subjected to subcritical loading, i.e.such that the elastic strain energy is not sufficient to cause instant failure, delay their failure because ofthe increase of available energy fracture created by the fluid drainage (Wang and Hong, 2012b).In this work, our aim is to shed light on the rate-dependent fracture process in the brain tissuecaused by fluid draining. Firstly, we present the results of fracture tests performed on porcine brainsamples. To this aim we use the wire cutting protocol, which is a well established method to measurethe fracture properties of soft materials, including viscous foodstuff and gels (Goh et al., 2005; Baldiet al., 2012; Forte et al., 2015). A computational model is then developed in the framework of finitestrain continuum mechanics, representing the large strain behaviour and fluid interaction through aporo-hyperelastic model (Simon, 1992). The numerical analyses are focused on the process of crackpropagation, which in the case of wire cutting develops after the initial stages of indentation and tissuerupture (Terzano, 2020). Through a simplified model of the fracture process in dissipative materials(Zhao, 2014), we are able to consider the energy dissipated by fluid-structure interaction as a functionof the loading rate. Finally, we provide a comparison is provided with the poroelastic behaviour of abiomimicking gel that was previously characterised by Forte et al. (2015).
2. Materials and methods
When measuring the fracture toughness of soft materials, traditional techniques based on stressintensity factors cannot be employed, since failure occurs when a large portion of the material is wellbeyond the limit of small strain elasticity. Toughness is hereby defined as the total amount of energyabsorption during deformation until fracture occurs (Huang et al., 2019). Wire cutting tests were4referred with respect to other available methods (such as, for instance, edge-notched tensile tests(Long and Hui, 2016)) because of the issues related to the extreme softness of the brain tissue, theeffect of self-weight and the impossibility of realising proper clamping. Porcine brain tissue sampleswere prepared, removing the cerebellum and separating the two hemispheres; each hemisphere was thenpositioned in the sample container with the frontal lobe facing upwards. The specimen would slowlyshift under gravity and approximately occupy a square prism of length 30 mm, width 30 mm and height50mm. Steel wires of diameter d w = 0 . , . , . , . v = 5mm s − , and the test with d w = 0 . v = 0 . − and v = 50mm s − .All tests were performed with a Biomomentum Mach-1 ™ mechanical testing system using a 1.5 N single-axis load cell, in a conditioned room at 19 ° C temperature (Forte et al., 2016). A schematic of theexperimental setup is shown in Fig. 1. load cellwirebrain sampleslit
Figure 1.
Wire cutting testing schematic. 1-col figure
A model is proposed to account for the rate-dependence observed during wire cutting tests onthe porcine brain samples, which can be extended to similar materials with a soft and wet porousmicrostructure. It is based on the following assumptions: (i) brain tissue and the biomimicking gels5re modelled as poro-hyperelastic materials; (ii) rate-dependent failure is described with a model ofthe fracture process based on the spatial separation of dissipative length scales and the definition ofa cohesive process zone; (iii) fracture in cutting is assimilated to the propagation of a far-field loadedcrack, depending on a geometric parameter (in this specific case, the wire diameter).
The brain tissue and the biomimicking gels are considered as biphasic materials, where a soft solidskeleton is saturated by biological fluids (Franceschini et al., 2006; Forte et al., 2017; Comellas et al.,2020). In this section we describe a poroelastic model at large strains, with specific focus on the equationsneeded for its numerical implementation in a finite element (FE) code. The theory of finite deformationcontinuum mechanics, as presented in standard textbooks on the subject, e.g. (Holzapfel, 2000), is thebackground in which the model is developed. For the sake of consistency with an updated Lagrangianframework in which the incremental solution strategy is implemented in the commercial FE code, thefield equations are referred to the current configuration. Accordingly, a material point of the biphasicmedium is identified by the position vector x ( X , t ), u S ≡ u = x − X ( x , t ) is the displacement of thispoint in the porous solid phase and u F defines the corresponding quantity for the pore fluid (Simon,1992) (Fig.2a). We also recall the decomposition of the spatial velocity gradient l = ∇ ˙ u = d + w , where d = sym( ∇ ˙ u ) is the symmetric rate of deformation tensor, w = − w T is the anti-symmetric spin tensorand ˙ u is the velocity of the solid phase. In a biphasic material, each phase in the current configuration is defined by a volume fraction n α = d v α / d v , where α = S, F corresponds, respectively, to the solid skeleton and the pore fluid.Assuming conditions of saturation, we establish the fundamental relationship n F + n S = 1 (Cheng,2016). In the following, we denote n = n F the porosity of the medium, which is correlated to thecurrent void ratio through e = n/ (1 − n ). The continuity mass equation for phase α reads (Ehlers and Throughout this section, ∇ ( • ) denotes the spatial gradient while ∇ · ( • ) is used for the spatial divergence operator.Italic is used for scalars, bold italic for vectors and bold roman for tensors. t ( n α ρ α ) + n α ρ α ∇ · ˙ u α = 0 (1)where D( • ) / D t is used for the material time derivative and ρ α is the effective density of each phase. Inthe solid skeleton, Eq.(1) provides the following relationship(1 − n )(1 − n ) = J S J − (2)where n is the porosity in the initial configuration, J is the volume ratio of the biphasic material and J S is instead referred to the solid skeleton. Notice that if the matrix is assumed incompressible, wehave J S = 1; however, this does not lead to macroscopic incompressibility, because volume change canoccur through changes in the volume fractions. Considering an incompressible fluid phase, Eq.(1) canbe rewritten as DD t n + n ∇ · ˙ u F = 0 (3)The strong form of the linear momentum balance for the biphasic material in quasi-static conditionsis provided by (Simon, 1992) ∇ · σ + ρ b = (4)where σ is the Cauchy stress tensor, ρ = (1 − n ) ρ S + nρ F is the homogenised density and b is the vectorof body forces per unit mass. The corresponding weak form can be written as (cid:90) Ω σ : δ e d v − (cid:90) Ω ρ b · δ u d v − (cid:90) ∂ Ω t ¯ t · δ u d a = 0 ∀ δ u (5)where δ e = sym ∇ δ u , with δ u being the virtual solid displacement, and ¯ t is the prescribed tractionvector on the boundary ∂ Ω t . In order to obtain the constitutive equations of the biphasic material in asuitable formulation for incremental Newton’s type numerical methods, a rate form should be introduced(Holzapfel, 2000). The first term in Eq.(5) represents the internal virtual work and taking its material7ime derivative we have ˙ δW int = (cid:90) Ω ( σ : ˙ δ e + ˙ σ : δ e ) d v (6)where ˙ δ e = − sym( ∇ δ u l ) (Holzapfel, 2000). In order to express the material rate of the Cauchy stress,we refer to the well-known concept of the effective stress in biphasic materials σ (cid:48) = σ + p F I (Biot, 1941),where p F = − / σ is the scalar pore pressure. Then, the second term of the parenthesis in Eq.(6)can be expressed by ˙ σ = σ (cid:48) J + w · σ + σ · w T − ˙ p F I (7)where we have introduced the objective Jaumann rate of the effective Cauchy stress σ (cid:48) J . The use ofthis specific objective rate is motivated by the implementation of the model in the commercial soft-ware Abaqus (Dassault Syst`emes SIMULIA, 2017). For completeness, we recall that a discretised andlinearised form of the previous equations is then needed, where the virtual displacement field is ap-proximated with suitable interpolation functions and variations are computed with respect to the fieldvariables of the problem —which in our case are represented by nodal displacements u and pore pressurevalues p F .Finally, we introduce the constitutive assumptions for the biphasic medium. The fluid flow throughthe porous skeleton is characterized by Darcy’s law, with an isotropic permeability tensor which re-mains unchanged during the deformation. In quasi-static conditions and neglecting inertia, Darcy’s lawcorrelates the rate of fluid volume to the pressure gradient˙ w = − κη F ∇ p F , (8)where ˙ w = n ( ˙ u F − ˙ u ) is the seepage velocity, representing the rate of fluid volume flowing through aunit normal area, κ is the intrinsic permeability and η F is the fluid viscosity (Cheng, 2016).The behaviour of the solid skeleton is specified by a hyperelastic isotropic strain energy function.Several studies related to brain mechanics (Budday et al., 2017a; Forte et al., 2017) have shown that8 modified one-term Ogden model provides optimal fit to experimental data (Ogden, 1972). The com-pressibility of the solid skeleton is implemented through the usual decomposition of the solid deformationinto isochoric and volumetric parts, such that the strain energy density is provided byΨ(¯ λ i ) + U ( J ) = 2 µα (cid:0) ¯ λ α + ¯ λ α + ¯ λ α − (cid:1) + 1 D ( J − (9)where µ , α and D = 2 /K S are material parameters and ¯ λ i = J − / λ i are the modified principalstretches (Holzapfel, 2000). The effective Cauchy stress tensor is split into its deviatoric and volumetriccomponents σ (cid:48) = s (cid:48) + p (cid:48) I (Selvadurai and Suvorov, 2016), with p (cid:48) = ∂U/∂J and (Connolly et al., 2019) s (cid:48) = J − β i ( n i ⊗ n i ) (10)where n i are the principal spatial directions and the stress coefficients are expressed by β i = ¯ λ i ∂ Ψ /∂ ¯ λ i − / λ j ∂ Ψ /∂ ¯ λ j (the summation rule applies to repeated indices). The last step required for the numericalimplementation in the updated Lagrangian framework is to make explicit the objective rate introducedin Eq.(7) through a spatial fourth-order elasticity tensor, such that σ (cid:48) J = c (cid:48) J : d (11)where c (cid:48) J is the spatial elasticity tensor defined in terms of the Jaumann rate of the Cauchy stress(Crisfield, 1997). The explicit formulation for a compressible hyperelastic model in terms of the principalstretches can be found, for instance, in the recent work by Connolly et al. (2019). The flow of interstitial fluid in the pores of the soft solid skeleton results in time-dependent deforma-tion and draining of the biphasic medium, according to a relaxation time which depends on the materialproperties (namely, the permeability) and the length of macroscopic observation (Hu and Suo, 2012).9 luid flow A - process zone B - crack-tip dissipative zone C - far-field region A B C l p l p l d cohesive zone σ max (b)(a) time t =0 time t xX u Ω Ω X ,x y y X ,x X ,x P O P d v S d v F d v=J d V b ∂ Ω t Figure 2. (a) Reference and current configurations of a biphasic continuum body. (b) Illustrative sketch ofthe fracture process in a poroelastic soft material. Shown in figure are (A) the process zone with radius l p , (B)the crack-tip dissipative zone with radius l p and (C) the far-field region. The enlarged view shows the processzone schematised with a rate-independent cohesive zone model. 2-col figure The analysis of rate-dependent fracture requires that poroelastic relaxation is considered as a source ofenergy dissipation correlated to crack propagation (Creton and Ciccotti, 2016), in which the length ofobservation is put in relation with some characteristic size of the fracture process.The model here proposed is based on the ideal situation illustrated in Fig. 2b, where a propagatingcrack in a semi-infinite body is shown with three regions where dissipative phenomena possibly occur(Long and Hui, 2016). Firstly, we consider damage phenomena occurring at the molecular scale, whichare condensed within the so-called process zone and account for the intrinsic toughness of the material.Within this region, energy dissipation might be affected by the loading rate and result in a tougheningeffect which is proportional to the velocity of crack propagation. For instance, Forte et al. (2015)explained the rate-dependent toughening in gelatins through a mechanism of fluid draining in the poreswithin the process zone. At a larger scale, dissipative terms are originated from relaxation in the bulkmaterial but become relevant to crack propagation only if they affect the crack-tip region, which webroadly define as the material affected by the vicinity of the crack. Usually, their effect is to preventthe crack driving force provided by external loading from being fully delivered to the crack (Qi et al.,10018). Finally, bulk processes in the far-field zone, which can cause macroscopic relaxation, are neglectedas they do not contribute to the fracture process. As a consequence of the proposed decomposition,we assume to split the fracture energy in two terms: the intrinsic term Γ o originating from the processzone, and an additional term Γ d due to energy dissipation in the crack-tip region affected by propagation(Zhao, 2014), so that we have Γ = Γ o + Γ d (12)Unlike the intrinsic toughness, Γ d cannot be treated as a property of the material because it isaffected by rate. With respect to the size of the crack-tip dissipative zone l d , we introduce a furtherhypothesis. In elastic soft materials, crack propagation is coupled with large deformations, and this ismotivated by the fact that the energy cost of creating new surfaces is comparable to the elastic strainenergy in the material. A length scale can be defined, known as elasto-adhesive length or length offlaw sensitive failure (Creton and Ciccotti, 2016; Chen et al., 2017), which separates by several orderof magnitudes most stiff solids from soft tissues. To a first approximation, it also represents the radiusof a blunted propagating crack in an elastic material (Creton and Ciccotti, 2016). We wish to clarifythat this has nothing to do with energy dissipation but simply characterises the concept of softness bya fracture mechanics point of view. Here we assume that the size of the crack-tip dissipative zone iscoincident with that of the large strain region, such that we have l d ≡ (cid:37) o = Γ o /E, l d (cid:29) l p (13)where E is the initial Young’s modulus of the material. In this work, it is assumed that energy dissipationis due to the drainage of fluid in the crack-tip region, that is, we are not considering viscoelasticity, whilethe intrinsic toughness Γ o is considered rate-independent. Conceptually, this is equivalent to employa cohesive process zone which enriches the continuum poro-hyperelastic model with a prescribed rate-independent stress-displacement relationship on the separation interface (Schwalbe et al., 2012) (see the11nlarged view in Fig. 2b). In line with our assumptions, the characteristic length l p of the cohesiveregion is much smaller than the size of the crack-tip dissipative region l d . Notice that, although severaltime-dependent cohesive models have been proposed, e.g. Noselli et al. (2016), our approach is toconsider that relaxation occurs outside the cohesive zone. Cutting involves deformation, friction and fracture. Here we focus on the steady-state, which isdeveloped after the initial stage of contact and indentation when the external work is converted intofracture energy for crack propagation (Terzano et al., 2018). Our aim is to establish the limits underwhich crack propagation in cutting can be compared to propagating a crack in symmetric far-fieldloading conditions. Furthermore, the plain strain assumption is introduced. A schematic of the modelis shown in Fig. 4a.With respect to crack propagation under remote loading, the finite size of the cutting tool adds anadditional length to the fracture process in cutting. It is speculated that the tool exerts some sort ofconstraint on the elastic blunting of the crack, which can be limited by the fact that the crack openingdisplacement is determined by the tool geometry (Hui et al., 2003; Zhang et al., 2019). In an elasticmaterial, the tip radius of a blunted crack (cid:37) o in condition of propagation is a material property andrepresents a characteristic length of the fracture process. It can be compared with the wire diameter d w in order to distinguish two different scenarios (Terzano, 2020): • for d w ≥ (cid:37) o , crack propagation happens as an autonomous process under symmetric mode-Iconditions, with a certain distance between the wire and the crack tip. The crack tip radius isdetermined by its natural value (cid:37) o = Γ o /E ; • for d w < (cid:37) o , the shape of the blunted crack is constrained by the wire, which touches the crack-tip. In this situation, the mechanism of propagation is different from that under remote loads andrequires a further input of external energy. In other terms, crack propagation is energy limited.12he analyses of fracture described in this work applies to wire cutting only when the first condition ismet. With this assumption, we have neglected the role of friction, which is known to affect the fracturetoughness of a material (Spagnoli et al., 2019). Furthermore, we have also considered that the criterionderived in an elastic situation is extended to the poro-hyperelastic material.
3. Results
Force-displacement curves obtained from wire cutting tests on the porcine brain tissue are illustratedin Fig.3a, for a wire diameter d w = 0 . o . To do so, we need toremove the contribution due to energy dissipation; typically, this means performing a fracture test atvery low loading rates, so that quasi-static conditions are assumed (Persson et al., 2005). The resultsare here elaborated according to the model proposed by Kamyab et al. (1998). Briefly, the steady-state cutting force F ss results from the force needed to open the crack and a contribution due to theformation of a flow zone around the bottom half of the wire, as shown in Fig.(3c). Friction produces acircumferential stress in this region but is neglected anywhere else. Then, the force per unit thicknessis proportional to the wire diameter, according to F ss t = Γ + (1 + f ) σ max d w (14)13 =0.5mm s -1 v =5mm s -1 v =50mm s -1 d w =0.16mm C u tt i ng f o r c e F / t [ N / m ] C u tt i ng f o r c e F / t [ N / m ] D [mm] Displacement D [mm] (a) (b) v =0.5mm s -1 v =5mm s -1 v =50mm s -1 Wire velocity v [mm/min] (log)0 0.2 0.4 0.6Wire diameter d w [mm] T oughne ss Г [ J m - ] ( l og ) [ N / m ] (d)(c) (e) C u tt i ng f o r c e F ss /t d w =0.16mmwhite matter failuregrey matter failurecrack propagation F max f max d w crackflow zone t D Г Г o
30 300 3000
Figure 3. (a) Cutting force versus displacement for d w = 0 . F ss /t as a function of the wire diameter. Thecontinuous linear fit is obtained for v = 5mms − . (e) Logarithmic plot of the intrinsic toughness as a functionof the insertion speed. 2-col figure where σ max has to be intended as a characteristic cohesive stress of the material, f is the frictionalcoefficient and t is the out-of-plane thickness of the sample.The steady-state force F ss /t obtained from the cutting experiments at v = 5mm s − is plotted asa function of the wire diameter in Fig.3d. Since a steady value cannot be easily recognised, F ss is14omputed as the average force corresponding to the onset of crack propagation observed in the tests.A linear fit is employed to extrapolate the force to zero diameter, such that, according to Eq.(14), thevalue of the fracture toughness is obtained. However, the experimental value of the force F ss /t for v = 0 . − is lower, suggesting that there might be an extra contribution due to energy dissipationresulting in Γ > Γ o . Lacking complete data for lower velocities, due to the complexity of realising propertests on the super-soft brain tissue, we then hypothesize that the same force-diameter slope applies toany insertion speed and extrapolate to d w = 0 the corresponding steady-state forces. These are shownon a logarithmic plane in Fig.(3e) and fitted with a linear interpolating function. The intercept withthe vertical axis, corresponding to a quasi-static value of the insertion speed, should provide the correctvalue of Γ o . Furthermore, the increase of toughness with speed follows a power-law, with exponentapproximately equal to 0.2. Due to the uncertainty in experimental data, we might assume that a valueof Γ o comprised between 0 . −
1J m − is a reasonable approximation. Indeed, this is the same order ofmagnitude of the toughness of biomimicking gelatins computed from wire cutting tests by Forte et al.(2015). Finite element analyses are employed to understand the origin of the rate-dependent fracture prop-erties observed in experiments. In order to reduce fracture in cutting to a problem of crack propagation,we first need to verify the hypothesis presented in Sect.2.2.3. We have modelled the steady state phaseas the insertion of a rigid circular wire into an edge-cracked body of width w and height 2 h , with initialcrack length c (Fig.4a). Since the wire extension in the out-of-plane direction is much larger than thethickness t of the samples, a plain strain assumption can be made. Due to symmetry, only half specimenis modelled with pertinent constraints imposed to the lower edge of the body; eight-node plane strainelements are employed, with a suitable refinement around the crack tip, which is artificially bluntedby taking an initial small radius of curvature. From analyses of mesh convergence, the smallest ele-15 ire y y y y w h c U U d w F ⇒ 𝜚=𝜚 o (a) (b) (c) J-integral contour
Wire diameter d w (mm) C r i t i c a l t i p r ad i u s 𝜚 ( mm ) d w =0.75 mm d w =0.25 mm d w =0.50 mm d w =0.125 mm y / c y / c Edge-crack -0.1 -0.05 000.040.080.120.16 𝜚 o Figure 4. (a) Sketch of the plane strain geometry employed to investigate the fracture process in wire cutting.Under the assumption of d w ≥ (cid:37) o , the equivalent model of an edge crack of length c , subjected to an openingremote displacement U , is obtained. The local system y , y is defined with respect to the moving crack tip. (b)Deformed crack profiles in the elastic brain tissue when J = Γ o , for an edge-crack (dashed line) and differentwire diameters d w . (c) Corresponding crack tip radius (cid:37) versus the wire diameter. The horizontal line is thenatural crack tip radius (cid:37) o . 2-col figure ment in the crack tip region is equal to 10 − h . The sample material is purely elastic, described bythe strain energy provided in Eq.(9). In such a case, the crack driving energy is correctly provided bythe J -integral, Eq.(16), such that the onset of crack propagation occurs when J = Γ o . The materialparameters implemented in the FE model are summarised in Table 1. The analyses were run with thequasi-static implicit solver of the commercial software Abaqus.Firstly, the case of an edge-crack subjected to far-field loading, by means of applied displacements U in the direction perpendicular to the crack, is considered. Then, we have studied the insertion of wireswith diameter d w = 0 . − y , y , the radius of the blunted crack can be expressed through the16 able 1. Mechanical parameters of the poro-hyperelastic model
Ogden’s parameters Brain and CH Gelatine (10% w/w) µ (Pa) 0 . · . · α -4.4 2.64 D Pa − . · − · − Hydraulic conductivity k m s − . · − . · − Fluid specific weight γ F (kN m − ) 9741 9741Initial void ratio e (%) 20 90Intrinsic toughness Γ o J m − Forte et al. (2017) Forte et al. (2015)radius (cid:37) of a circle fitting the profile within a distance equal to 10 − c from the crack tip. Plots of thedeformed crack when J = Γ o are shown in Fig.4b, suggesting that the crack blunting with wires ofvarious diameters is almost equivalent to the edge-crack subjected to remote loading. The transition toconstrained blunting seems to occur when d w = 0 . (cid:37) o = 2Γ o /E , which forthe soft tissues of our study is in the order of 1 . · − m. By plotting the critical crack tip radius (cid:37) against the wire diameter d w Fig.(4c), we notice that it is approximately equal to (cid:37) o when d w ≥ (cid:37) o .Below this limit, we hypothesise that the tip radius scales with the wire diameter (hence the slope 1/2shown in the plots). In conclusion, we can assume that, in the materials under consideration, steadystate cutting is equivalent to a problem of crack propagation when the wire diameter is d w ≥ . Retaining the assumption of autonomous crack propagation, we can study the rate-dependent frac-ture in an equivalent edge-cracked model with applied remote displacements. The geometry is illustratedin Fig.5a: it consists of a large rectangular sample of height 2 h = 50mm and width w = 20mm, con-taining an edge-crack of length c = 1mm. Normal displacements are applied to the top and bottomboundaries such that the strain rate is constant, that is U = [exp( ˙ εt ) − h , where ε = ln[( h + U ) /h ] is17 -1 =10 min -1 =100 min -1 crack crack tip zone p F = λ =1+ U / hλ =1+ U / h λ = λ o hhw pore pressure p F p F = p F = c (a) (b) brain gelatine brain gelatine Strain rate (min -1 , log) 1 10 100 10000.7511.251.5 C r i t i c a l s t r e t c h -1 , log) λ / λ o (c) (d)(a) gelatinebrain tissue l d l d l d Figure 5. (a) Model of the edge-cracked sample for the analyses of poroelastic fracture. Blue edges are thosewith drained boundary conditions. (b) Contours of the pore pressure p F at constant stretch λ = λ o for differentstrain rates, in the brain tissue and the gelatine studied by Forte et al. (2015).(c) Strain energy per unit areain the crack-tip region, normalised by the fracture toughness Γ o , and (d) applied stretch, normalised by thecritical stretch in quasi-static conditions λ o . Both are plotted as a function of the strain rate ˙ ε (logarithmicplot). 2-col figure the true strain in the direction normal to the crack line. The stretch ratio is defined by λ = 1 + U/h .Two materials, the brain tissue and the gelatine studied by Forte et al. (2015), are described with theporo-hyperelastic model presented in Sect.2.2. The properties are summarised in Table 1. Notice thatthe hydraulic conductivity k is employed in place of the permeability κ , to which is related by means of18 = κγ F /η F . The finite element mesh is built with four-node quadrilateral plane strain hybrid elementswith additional degrees of freedom for the pore pressure p F . Boundary conditions are specified in termsof displacements (top and bottom forces are prevented from lateral motion), and in addition on thepore pressure degree of freedom. A condition of draining, enforced by setting the pore pressure equalto zero, is specified for the vertical free edges and the edge-crack surfaces in contact with atmosphericpressure (Fig.5a). The reference porosity n needs to be specified as initial condition through the voidratio e . The analyses were run with the implicit solver of the commercial software Abaqus. A transientfluid-stress diffusion analysis is required to simulate fluid flow through the porous material, where theaccuracy of the solution is governed by the maximum pore pressure change allowed in an increment.Different values have been considered for the best compromise between accuracy and efficiency.The main purpose of the analyses is to understand how fluid draining in the crack-tip region affectsthe onset of crack propagation. In other terms, we are considering the effect of dissipation and of theloading rate on the crack driving energy, whereas the fracture toughness is assumed equal to Γ o . Thecritical condition is then defined by J ( ˙ ε ) = Γ o (15)where J denotes the J -integral, which in fracture mechanics is used to compute the driving force for crackgrowth. Due to coupling between deformation of the solid network and fluid diffusion, the J -integral ispath-dependent since it includes poroelastic dissipation, unless a vanishing small contour surroundingthe crack tip is considered (Wang and Hong, 2012b). In this work we compute the J -integral in thebiphasic material according to (Shih et al., 1986) J = (cid:90) C (cid:18) Ψ n − ∂u i ∂x σ (cid:48) ij n j (cid:19) d s, (16)where Ψ is the strain energy density and n j is the unit vector normal to a contour C enclosing thecrack tip (Fig.4a). The results presented below are obtained considering a contour that surrounds the19rack-tip region of radius l d .A preliminary analysis on the elastic material is employed to investigate the quasi-static situation. Insuch a case, the critical condition evaluated through Eqs.(15)-(16) provides the stretch λ o , correspondingto the onset of crack propagation. Results in the poroelastic materials are illustrated in Fig.5b, where weshow the contours of the fluid pressure p F , for three different strain rates ˙ ε , when λ = λ o . The enlargedregion is the crack-tip dissipative zone, whose radius, from Eq.(13), is approximated to l d ∼ − m. Thered areas correspond to the drained or relaxed condition ( p F = 0) whereas the blue regions are affectedby fluid flowing in the pores. It can be seen that, independently from the rate, the greater permeabilityof gelatins allows for a rapid draining of the whole crack-tip region. On the contrary, it appears that fluidtakes a longer time to drain the same area in the brain tissue, where permeability is much lower. Sincefluid draining is a dissipative process, it is reasonable to assume that crack propagation is affected by thephenomenon, at least in the brain tissue. Keeping in mind the limitations in the use of the J -integral, inFig.5c we present the normalised energy at constant stretch λ o for different strain rates. The observedbehaviour can be better comprehended by plotting the normalised stretch when J = Γ o , Fig.5d. Asexpected, no difference with respect to the elastic quasi-static situation emerges in the gelatine, whichtherefore behaves as an elastic relaxed material. The situation looks different in the brain tissue, whereboth the strain energy and the critical stretch are affected by rate. Notice that we cannot consider thesestretches as the real ultimate stretches of the material; nevertheless, the results shown in Fig.5d suggesta toughening effect due to fluid draining in the brain tissue.
4. Discussion
The phenomenon of poroelastic relaxation is characterised by a long-range motion of interstitial fluidin the solid network (Wang and Hong, 2012a). The characteristic time of relaxation depends on thelength of observation, and this implies that samples with macroscopic sizes, such as those employed inthe experiments, require relatively long times compared, for instance, to viscoelastic relaxation (Hu andSuo, 2012). However, with respect to the presence of a crack the situation is completely different. The20opic of fracture in the brain tissue is not well documented; for this reason, we have developed our modelbased on the observations in hydrogels, which can be considered as benchmark examples of soft porousmaterials.In the crack tip zone of soft porous tissues, high uniaxial tensile stresses trigger a take-up of fluidfrom the far-field region (Long and Hui, 2015). Since these high stress gradients are confined to a smallregion, whose extension in the soft tissues under investigation is in the order of 10 − m, fluid drainingof the crack tip zone is a relatively quick phenomenon. Evidently, it might happen that viscoelasticrelaxation occurs at the same time, but this aspect has not been considered in the present work. Thekey observation is that material relaxation is relevant with respect to rate-dependent fracture in relationto the crack-tip dissipative region. In our analysis, we have explained this fact with separation of lengthscales, through which we can neglect the dissipative phenomena in the bulk; furthermore, we haveassumed that the rate-independent toughness threshold is originated within the process zone, whoseextension is in the scale of nanometers. There are models that explained fluid-related toughening basedon the fluid draining in the process zone (Forte et al., 2015). Interestingly, the model by Forte et al.(2015) was proposed to explain the toughening of gelatins, which cannot be captured by the analysespresented in this work (Figs.5b-d).With specific consideration of fluid-related effects, the mechanism of crack-tip draining is illustratedin Fig.6. Initially the solid is saturated and there is only an infinitesimal zone close to the tip of theexisting crack where the fluid pressure is zero; as time goes by, this region increases in extent until thewhole crack-tip region is drained. According to Biot’s theory, the pressure-driven fluid flow is a diffusiveprocess. The time of poroelastic relaxation can be defined as the time needed to drain an area of radius l d close to the crack-tip (Hui et al., 2013) t d = l d /D F (17)where D F is the diffusion coefficient, which depends on the permeability, the fluid viscosity and the elasticproperties of the solid. To a first approximation, in linear poroelasticity and plane strain conditions we21 rained undrainedcrack tipcrack tipfast rateslow rate l d l d Figure 6.
Sketch of the crack-tip dissipative region and influence of fluid draining at fast and slow strain rates˙ ε . 1-col figure have D F = 2 µ (1 − ν ) κ/ (1 − ν ) η F , where µ and ν are, respectively, the shear modulus and Poisson’sratio of the solid material (Hui et al., 2013). Ideally, we can distinguish two limit situations: • fast strain rates or reduced permeability: fluid diffusion is too slow to be effective and the cracktip region is saturated at the instant corresponding to the onset of propagation. With respect tofracture, the material behaves as a soft incompressible solid; • slow strain rates or high permeability: the draining process is fast, therefore a large drained regionsurrounds the crack-tip zone. Here we observe an extensive process of draining but this becomesineffective with respect to fracture. The material behaves as a soft compressible solid.Intermediate situations are those in which fluid drains the crack-tip region in a time range comparableto that leading to crack propagation, hence causing a dissipative phenomenon that may produce anenhancement of the material toughness. This effect was correctly described by our model in the braintissue. However, there are some limits in the procedure adopted, which are discussed below, that leavespace for further work on this topic.In the procedure through which we have derived the toughness of brain tissue from the experimental22ire cutting force, we have assumed that the slope of the steady state cutting force F SS versus the wirediameter d w is the same for different velocities (Fig.3d). However, wire cutting analyses on biomim-icking gelatins by Forte et al. (2015) have shown that the slope of the interpolating function increaseswith the insertion velocity, although it tends to become constant at lower velocities. Furthermore, theextrapolation procedure to infer the intrinsic toughness Γ o from Eq.(14) was based on the assumptionthat v = 1 mm min − is a reasonable speed for quasi-static conditions. To a first approximation, wemight relate the quasi-static threshold to the process of fluid diffusion, which in turn depends on thepermeability of the material. The employed value was derived from similar observations by Forte et al.(2015) on gelatins. However, considering the large difference with the permeability of brain tissue, ourassumption needs to be verified against further experimental observations. This point brings us to acentral aspect in modelling fluid-related effects in soft tissues: the issue of accurately measuring andmodelling permeability. For the sake of simplicity, we have adopted the hypothesis of material isotropy:however, while it seems to be a valid assumption for the elasticity of the brain tissue (Budday et al.,2020), diffusion or permeability properties are remarkably anisotropic, in particular in white matterregions characterised by axonal structures. In addition, brain tissue permeability can be modified sub-stantially under loading by swelling and additional coupling with the local tissue deformation (Jamalet al., 2020).The model that we have developed isolated the stage of crack propagation, leaving aside the wholeprocess of contact and indentation that occurs in wire cutting. Ideally, we could have simulated thecomplete cutting process directly through the finite element model and use cohesive elements to simulatethe process of propagation. This approach has successfully modelled needle penetration in soft elasticmaterials (Oldfield et al., 2013; Terzano et al., 2020) and rate-dependent wire cutting of viscous food(Goh et al., 2005; Skamniotis and Charalambides, 2020). However, to calibrate the cohesive model forthe brain tissue would require the characterisation of its frictional behaviour (Casanova et al., 2014),which is known to affect the fracture toughness of the material (Duncan et al., 2020). In addition,the role of fluid on the frictional contact between wire and tissue should also be considered (Reale and23unn, 2017), possibly implying some effect of lubrication which, at the present time, is unknown. TheFE model provided meaningful results on poroelastic toughening with rate (Figs.5c-d) but we haveno means to establish a direct confrontation with experimental data. Experiments revealed a rate-dependent toughening in terms of the velocity of wire insertion (Figs.3e), which in the steady state canbe reasonably considered to coincide with the crack velocity. In the numerical model we have insteadexplored the effect of the strain rate on the onset of crack propagation, but we cannot establish ananalytical relationship between the strain rate and the crack propagation velocity. It comes naturallyto think that higher strain rates result in faster crack propagation, although this might hold only belowa certain limit, as shown for instance in fracture tests on hydrogels (Mayumi et al., 2016).Finally, although our poro-hyperelastic model is able to couple large deformation and fluid flow,more advanced models considering the full coupling between the solvent diffusion and tissue swellingmight be required, e.g. Hong et al. (2008); Bouklas et al. (2015a); Chen et al. (2020); Brighenti andCosma (2020). This would also allow us to implement a modified definition of the J -integral proposedfor swelling materials, which is path-independent and computes the transient energy energy release rateby separating the energy lost in diffusion from the energy available to drive crack growth (Yang et al.,2006; Bouklas et al., 2015b).
5. Conclusion
Testing the fracture properties of super soft tissues through standard tensile specimens is a complextask. For this reason, wire cutting was here employed to analyse the influence of rate on the fractureenergy of brain tissue. The experimental data show an evident increase of the cutting force with therate of insertion, suggesting that some form of energy dissipation affects the cutting process. In thiswork, we speculate that the rate-dependent toughening is due to poroelastic dissipation in the vicinity ofthe crack that is propagated ahead of the wire. We have proposed a numerical model which consideredthe brain tissue as a biphasic material. Through finite elements analyses of an edge-cracked sample,subjected to remote loading with varying strain rate, we have shown how the process of fluid draining24n the crack-tip region might affect the fracture toughness of the material. We can then summarise themain findings: • the analysis of wire cutting experimental data suggests a power-law increase of fracture toughnesswith the rate of insertion, as already observed in biomimicking gelatins; • we have identified a length scale which distinguishes the fracture process of cutting from crackpropagation under remote loading. Specifically, below a certain wire diameter crack propagationbecomes unstable and the shape of the crack is constrained. Interestingly, this was observedexperimentally on hydrogels (Baldi et al., 2012) and was motivated by the reduced stiffness ofsuch materials. More correctly, we are able to say that it depends on the competition between thecost of creating new surfaces and the elastic strain energy of the material; • the finite element analyses of the fracture process in the poroelastic material have confirmed thetoughening effect with the rate of applied loading. According to our poro-hyperelastic model, sucha contribution is chiefly controlled by the value of the intrinsic permeability of the material.This work has purposely neglected the dissipative behaviour provided by viscoelasticity in order tofocus on fluid-related effects. Future work will be dedicated to extend the proposed model to coupledviscoelasticity and fluid diffusion. In the context of fracture, accurate models should specifically targetthe rate-sensitivity of the process zone. By a computational point of view, cohesive models might stillbe the ideal candidates to include energy dissipation through a time-dependent cohesive law. Our viewis that they should be developed on the ground of a micromechanical description of the disintegratingmaterial ahead of the crack tip. In particular, we envisage that further research is needed to characterisethe effect of water diffusion on mechanical deformation by a micromechanical point of view.25 . Acknowledgements A. E. Forte acknowledges the support received from the European Union’s Horizon 2020 research andinnovation programme under the Marie Sk(cid:32)lodowska-Curie grant agreement No 798244. The authors alsoacknowledge the financial support from EDEN2020 project funded by the European Union’s Horizon2020 research and innovation programme under grant agreement No 688279. D. Dini would like toacknowledge the support received from the UKRI Engineering and Physical Sciences Research Council(EPSRC) via his Established Career Fellowship EP/N025954/1.
References
Ahmadzadeh, H., Smith, D.H., Shenoy, V.B., 2014. Viscoelasticity of Tau Proteins Leads to Strain Rate-Dependent Breaking of Microtubules during Axonal Stretch Injury: Predictions from a MathematicalModel. Biophysical Journal 106, 1123–1133. doi: .Arroyo, M., Trepat, X., 2017. Hydraulic fracturing in cells and tissues: fracking meets cell biology.Current Opinion in Cell Biology 44, 1–6. doi: .Ashammakhi, N., Ahadian, S., Darabi, M.A., El Tahchi, M., Lee, J., Suthiwanich, K., Sheikhi, A., Dok-meci, M.R., Oklu, R., Khademhosseini, A., 2019. Minimally Invasive and Regenerative Therapeutics.Advanced Materials 31, 1804041. doi: .Baldi, F., Bignotti, F., Peroni, I., Agnelli, S., Ricc`o, T., 2012. On the measurement of the fractureresistance of polyacrylamide hydrogels by wire cutting tests. Polymer Testing 31, 455–465. doi: .Baumberger, T., Caroli, C., Martina, D., 2006. Solvent control of crack dynamics in a reversible hydrogel.Nature Materials 5, 552–555. doi: , arXiv:0605395 .26iot, M.A., 1941. General theory of three-dimensional consolidation. Journal of Applied Physics 12,155–164. doi: .Bouklas, N., Landis, C.M., Huang, R., 2015a. A nonlinear, transient finite element method for coupledsolvent diffusion and large deformation of hydrogels. Journal of the Mechanics and Physics of Solids79, 21–43. doi: .Bouklas, N., Landis, C.M., Huang, R., 2015b. Effect of Solvent Diffusion on Crack-Tip Fields and DrivingForce for Fracture of Hydrogels. Journal of Applied Mechanics 82, 081007. doi: .Brighenti, R., Cosma, M.P., 2020. Swelling mechanism in smart polymers responsive to mechano-chemical stimuli. Journal of the Mechanics and Physics of Solids 143, 104011. doi: .Budday, S., Ovaert, T.C., Holzapfel, G.A., Steinmann, P., Kuhl, E., 2020. Fifty Shades of Brain: AReview on the Mechanical Testing and Modeling of Brain Tissue. Archives of Computational Methodsin Engineering 27, 1187–1230. doi: .Budday, S., Sommer, G., Birkl, C., Langkammer, C., Haybaeck, J., Kohnert, J., Bauer, M., Paulsen, F.,Steinmann, P., Kuhl, E., Holzapfel, G.A., 2017a. Mechanical characterization of human brain tissue.Acta Biomaterialia 48, 319–340. doi: .Budday, S., Sommer, G., Holzapfel, G.A., Steinmann, P., Kuhl, E., 2017b. Viscoelastic parameteridentification of human brain tissue. Journal of the Mechanical Behavior of Biomedical Materials 74,463–476. doi: .Casanova, F., Carney, P.R., Sarntinoranont, M., 2014. In vivo evaluation of needle force and frictionstress during insertion at varying insertion speed into the brain. Journal of Neuroscience Methods237, 79–89. doi: .27hen, C., Wang, Z., Suo, Z., 2017. Flaw sensitivity of highly stretchable materials. Extreme MechanicsLetters 10, 50–57. doi: .Chen, S., Huang, R., Ravi-Chandar, K., 2020. Linear and nonlinear poroelastic analysis of swelling anddrying behavior of gelatin-based hydrogels. International Journal of Solids and Structures 195, 43–56.doi: .Cheng, A.H.D., 2016. Poroelasticity. Theory and Applications of Transport in Porous Media, SpringerInternational Publishing, New York.Cloots, R.J.H., van Dommelen, J.A.W., Nyberg, T., Kleiven, S., Geers, M.G.D., 2011. Micromechanicsof diffuse axonal injury: influence of axonal orientation and anisotropy. Biomechanics and Modelingin Mechanobiology 10, 413–422. doi: .Comellas, E., Budday, S., Pelteret, J.P., Holzapfel, G.A., Steinmann, P., 2020. Modeling the porousand viscous responses of human brain tissue behavior. Computer Methods in Applied Mechanics andEngineering 369, 113128. doi: .Connolly, S.J., Mackenzie, D., Gorash, Y., 2019. Isotropic hyperelasticity in principal stretches: explicitelasticity tensors and numerical implementation. Computational Mechanics 64, 1273–1288. doi: .Creton, C., Ciccotti, M., 2016. Fracture and adhesion of soft materials: a review. Reports on Progressin Physics 79, 046601. doi: .Crisfield, M.A., 1997. Non-linear Finite Element Analysis of Solids and Structures: Advanced topics.Wiley, Chichester.Dassault Syst`emes SIMULIA, 2017. Abaqus 2017, Documentation.28uncan, T.T., Sarapas, J.M., Defante, A.P., Beers, K.L., Chan, E.P., 2020. Cutting to measure theelasticity and fracture of soft gels. Soft Matter doi: .Ehlers, W., Eipper, G., 1999. Finite Elastic Deformations in Liquid-Saturated and Empty PorousSolids, in: Porous Media: Theory and Experiments. Springer Netherlands, Dordrecht. volume 34, pp.179–191.Ehlers, W., Wagner, A., 2015. Multi-component modelling of human brain tissue: a contributionto the constitutive and computational description of deformation, flow and diffusion processes withapplication to the invasive drug-delivery problem. Computer Methods in Biomechanics and BiomedicalEngineering 18, 861–879. doi: .El Sayed, T., Mota, A., Fraternali, F., Ortiz, M., 2008. Biomechanics of traumatic brain injury. Com-puter Methods in Applied Mechanics and Engineering 197, 4692–4701. doi: .Forte, A.E., D’Amico, F., Charalambides, M.N., Dini, D., Williams, J.G., 2015. Modelling and experi-mental characterisation of the rate dependent fracture properties of gelatine gels. Food Hydrocolloids46, 180–190. doi: .Forte, A.E., Galvan, S., Dini, D., 2018. Models and tissue mimics for brain shift simulations. Biome-chanics and Modeling in Mechanobiology 17, 249–261. doi: .Forte, A.E., Galvan, S., Manieri, F., Rodriguez y Baena, F., Dini, D., 2016. A composite hydrogel forbrain tissue phantoms. Materials & Design 112, 227–238. doi: .Forte, A.E., Gentleman, S.M., Dini, D., 2017. On the characterization of the heterogeneous mechanicalresponse of human brain tissue. Biomechanics and Modeling in Mechanobiology 16, 907–920. doi: . 29ranceschini, G., Bigoni, D., Regitnig, P., Holzapfel, G.A., 2006. Brain tissue deforms similarly tofilled elastomers and follows consolidation theory. Journal of the Mechanics and Physics of Solids 54,2592–2620. doi: .Garc´ıa, J.J., Smith, J.H., 2009. A Biphasic Hyperelastic Model for the Analysis of Fluid andMass Transport in Brain Tissue. Annals of Biomedical Engineering 37, 375–386. doi: .Goh, S., Charalambides, M.N., Williams, J.G., 2005. On the mechanics of wire cutting of cheese.Engineering Fracture Mechanics 72, 931–946. doi: .Guimar˜aes, C.F., Gasperini, L., Marques, A.P., Reis, R.L., 2020. The stiffness of living tissuesand its implications for tissue engineering. Nature Reviews Materials 5, 351–370. doi: .Haldar, K., Pal, C., 2018. Rate dependent anisotropic constitutive modeling of brain tissue undergoinglarge deformation. Journal of the Mechanical Behavior of Biomedical Materials 81, 178–194. doi: .Holzapfel, G.A., 2000. Nonlinear Solid Mechanics: A Continuum Approach for Engineering Science.Wiley, Chichester.Hong, W., Zhao, X., Zhou, J., Suo, Z., 2008. A theory of coupled diffusion and large deformation inpolymeric gels. Journal of the Mechanics and Physics of Solids 56, 1779–1793. doi: .Hosseini-Farid, M., Ramzanpour, M., McLean, J., Ziejewski, M., Karami, G., 2020. A poro-hyper-viscoelastic rate-dependent constitutive modeling for the analysis of brain tissues. Journal of theMechanical Behavior of Biomedical Materials 102, 103475. doi: .30u, Y., Suo, Z., 2012. Viscoelasticity and poroelasticity in elastomeric gels. Acta Mechanica SolidaSinica 25, 441–458. doi: .Huang, W., Restrepo, D., Jung, J., Su, F.Y., Liu, Z., Ritchie, R.O., McKittrick, J., Zavattieri, P., Ki-sailus, D., 2019. Multiscale Toughening Mechanisms in Biological Materials and Bioinspired Designs.Advanced Materials 31, 1901561. doi: .Hui, C.Y., Jagota, A., Bennison, S.J., Londono, J.D., 2003. Crack blunting and the strength of softelastic solids. Proceedings of the Royal Society of London. Series A: Mathematical, Physical andEngineering Sciences 459, 1489–1516. doi: .Hui, C.Y., Long, R., Ning, J., 2013. Stress Relaxation Near the Tip of a Stationary Mode I Crack in aPoroelastic Solid. Journal of Applied Mechanics 80. doi: .Jamal, A., Mongelli, M.T., Vidotto, M., Madekurozwa, M., Bernardini, A., Overby, D.R., De Momi,E., Rodriguez y Baena, F., Sherwood, J.M., Dini, D., 2020. Infusion Mechanisms in Brain WhiteMatter and its Dependence of Microstructure: An Experimental Study of Hydraulic Permeability.IEEE Transactions on Biomedical Engineering , 1–1doi: .Jin, X., Zhu, F., Mao, H., Ming, S., Yang, K.H., 2013. A comprehensive experimental study on materialproperties of human brain tissue. Journal of Biomechanics 46, 2795–801. doi: .Kamyab, I., Chakrabarti, S., Williams, J.G., 1998. Cutting cheese with wire. Journal of MaterialsScience 33, 2763–2770. doi: .Long, R., Hui, C.Y., 2015. Crack tip fields in soft elastic solids subjected to large quasi-static deformation- A review. Extreme Mechanics Letters 4, 131–155. doi: .Long, R., Hui, C.Y., 2016. Fracture toughness of hydrogels: measurement and interpretation. SoftMatter 12, 8069–8086. doi: .31ak, A.F., 1986. Unconfined compression of hydrated viscoelastic tissues: A biphasic poroviscoelasticanalysis. Biorheology 23, 371–383. doi: .Mayumi, K., Guo, J., Narita, T., Hui, C.Y., Creton, C., 2016. Fracture of dual crosslink gels withpermanent and transient crosslinks. Extreme Mechanics Letters 6, 52–59. doi: .Mow, V.C., Kuei, S.C., Lai, W.M., Armstrong, C.G., 1980. Biphasic Creep and Stress Relaxation ofArticular Cartilage in Compression: Theory and Experiments. Journal of Biomechanical Engineering102, 73–84. doi: .Naassaoui, I., Ronsin, O., Baumberger, T., 2018. A poroelastic signature of the dry/wet state of a cracktip propagating steadily in a physical hydrogel. Extreme Mechanics Letters 22, 8–12. doi: .Nicholson, C., 2001. Diffusion and related transport mechanisms in brain tissue. Reports on Progressin Physics 64, 815–884. doi: .Noselli, G., Lucantonio, A., McMeeking, R.M., DeSimone, A., 2016. Poroelastic toughening in polymergels: A theoretical and numerical study. Journal of the Mechanics and Physics of Solids 94, 33–46.doi: , arXiv:1605.02000 .Ogden, R.W., 1972. Large Deformation Isotropic Elasticity - On the Correlation of Theory and Exper-iment for Incompressible Rubberlike Solids. Proceedings of the Royal Society of London. Series A:Mathematical, Physical and Engineering Sciences 326, 565–584. doi: .Oldfield, M.J., Dini, D., Giordano, G., Rodriguez y Baena, F., 2013. Detailed finite element modellingof deep needle insertions into a soft tissue phantom using a cohesive approach. Computer Methods inBiomechanics and Biomedical Engineering 16, 530–543. doi: .32ersson, B.N.J., Albohr, O., Heinrich, G., Ueba, H., 2005. Crack propagation in rubber-like materials.Journal of Physics: Condensed Matter 17, 1071–1142. doi: .Qi, Y., Caillard, J., Long, R., 2018. Fracture toughness of soft materials with rate-independent hysteresis.Journal of the Mechanics and Physics of Solids 118, 341–364. doi: .Qian, L., Zhao, H., Guo, Y., Li, Y., Zhou, M., Yang, L., Wang, Z., Sun, Y., 2018. Influence of strain rateon indentation response of porcine brain. Journal of the Mechanical Behavior of Biomedical Materials82, 210–217. doi: .Reale, E.R., Dunn, A.C., 2017. Poroelasticity-driven lubrication in hydrogel interfaces. Soft Matter 13,428–435. doi: .de Rooij, R., Kuhl, E., 2016. Constitutive Modeling of Brain Tissue: Current Perspectives. AppliedMechanics Reviews 68, 010801. doi: .Schwalbe, K.H., Scheider, I., Cornec, A., 2012. Guidelines for Applying Cohesive Models to the DamageBehaviour of Engineering Materials and Structures. Springer, Heidelberg.Selvadurai, A.P., Suvorov, A.P., 2016. Coupled hydro-mechanical effects in a poro-hyperelastic material.Journal of the Mechanics and Physics of Solids 91, 311–333. doi: .Shih, C.F., Moran, B., Nakamura, T., 1986. Energy release rate along a three-dimensional crack front ina thermally stressed body. International Journal of Fracture 30, 79–102. doi: .Simon, B.R., 1992. Multiphase Poroelastic Finite Element Models for Soft Tissue Structures. AppliedMechanics Reviews 45, 191–218. doi: .Skamniotis, C., Charalambides, M.N., 2020. Development of computational design tools for char-acterising and modelling cutting in ultra soft solids. Extreme Mechanics Letters 40, 100964.doi: . 33pagnoli, A., Brighenti, R., Terzano, M., Artoni, F., 2019. Cutting resistance of soft materials: Effectsof blade inclination and friction. Theoretical and Applied Fracture Mechanics 101, 200–206. doi: .Stastna, M., Tenti, G., Sivaloganathan, S., Drake, J.M., 1999. Brain biomechanics: consolidationtheory of hydrocephalus. Variable permeability and transient effects. Canadian Applied MathematicsQuarterly 7, 93–109.Terzano, M., 2020. Fracture processes in indentation and cutting of soft biomaterials. Phd thesis.University of Parma.Terzano, M., Dini, D., Rodriguez y Baena, F., Spagnoli, A., Oldfield, M., 2020. An adaptive finiteelement model for steerable needles. Biomechanics and Modeling in Mechanobiology doi: .Terzano, M., Spagnoli, A., St˚ahle, P., 2018. A fracture mechanics model to study indentation cutting.Fatigue & Fracture of Engineering Materials and Structures 41, 821–830. doi: .Wang, X., Hong, W., 2012a. A visco-poroelastic theory for polymeric gels. Proceedings of the RoyalSociety of London. Series A: Mathematical, Physical and Engineering Sciences 468, 3824–3841. doi: .Wang, X., Hong, W., 2012b. Delayed fracture in gels. Soft Matter 8, 8171. doi: .Yang, F., Wang, J., Chen, D., 2006. The energy release rate for hygrothermal coupling elastic materials.Acta Mechanica Sinica 22, 28–33. doi: .Zhan, W., Rodriguez y Baena, F., Dini, D., 2019. Effect of tissue permeability and drug diffusionanisotropy on convection-enhanced delivery. Drug Delivery 26, 773–781. doi: . 34hang, B., Shiang, C.S., Yang, S., Hutchens, S., 2019. Y-Shaped Cutting for the SystematicCharacterization of Cutting and Tearing. Experimental Mechanics 59, 517–529. doi: .Zhao, X., 2014. Multi-scale multi-mechanism design of tough hydrogels: building dissipation intostretchy networks. Soft Matter 10, 672–687. doi:10.1039/C3SM52272E