## Abstract

In some soft biological structures such as brain and fat tissues, strong experimental evidence suggests that the shear modulus increases significantly under increasing compressive strain, but not under tensile strain, whereas the apparent Young's elastic modulus increases or remains almost constant when compressive strain increases. These tissues also exhibit a predominantly isotropic, incompressible behaviour. Our aim is to capture these seemingly contradictory mechanical behaviours, both qualitatively and quantitatively, within the framework of finite elasticity, by modelling a soft tissue as a homogeneous, isotropic, incompressible, hyperelastic material and comparing our results with available experimental data. Our analysis reveals that the Fung and Gent models, which are typically used to model soft tissues, are inadequate for the modelling of brain or fat under combined stretch and shear, and so are the classical neo-Hookean and Mooney–Rivlin models used for elastomers. However, a subclass of Ogden hyperelastic models are found to be in excellent agreement with the experiments. Our findings provide explicit models suitable for integration in large-scale finite-element computations.

## 1. Introduction

Obtaining reliable constitutive models for the behaviour of tissues under loads is of the utmost importance when studying the response and evolution of organs in physiological and pathological conditions. For instance, the computational analysis of traumatic brain injury owing to shocks or blast waves in sports, combat or accidents relies on large finite-element codes based on the constitutive properties of brain tissues. Similarly, an understanding of how brain tumours change the mechanical and neurological environment during growth depends on the mechanical responses of both healthy tissue and tumours [1]. The response of adipose tissue to external loads is also a growing area of interest in clinical research, for example in treating patients with impaired mobility and in the pharmaceutical industry, particularly for the design of needle-free drug-delivery systems [2].

Recent experimental evidence [3–5] shows that soft biological tissues such as brain, gliomas, liver and fat have some unusual mechanical properties under loads, namely

(i) the shear modulus increases sharply as compression in the direction orthogonal to the shear direction increases;

(ii) the shear modulus remains almost constant or may decrease as tension in the direction orthogonal to the shear direction increases; and

(iii) the elastic modulus increases or remains almost constant when compression increases.

In particular, the shear modulus of normal brain can be increased nearly four times by compressive stresses. In addition, although at low strains, glioma brain tumours have similar elastic moduli to normal brain tissue (unlike other tumour types arising in breast tissue, for example [6]), at large strains, glioma tissue stiffens more strongly under compression than normal brain. However, while the shear modulus increases significantly when axial strain increases, the elastic modulus increases only slightly or not at all under increasing axial strain. During experimental tests, tissue samples exhibited a predominantly isotropic behaviour and their volume was reported to remain virtually constant.

Our aim is to capture these seemingly contradictory mechanical behaviours, both qualitatively (theoretically) and quantitatively (numerically), within the framework of finite elasticity, by modelling a soft tissue as a homogeneous, isotropic, incompressible, hyperelastic material. First, we demonstrate analytically that, in large strain deformations, conditions (i)–(iii) can be satisfied simultaneously by Mooney–Rivlin models, but not by the neo-Hookean, Fung and Gent models. The neo-Hookean model can be derived from first principles and is suitable for materials with entropic elasticity and a Gaussian distribution of chains with quadratic strain energy. While the neo-Hookean model can be seen as a general second-order approximation of a strain–energy density, the Mooney–Rivlin model for incompressible systems is its third-order approximation and is known to be better suited than the neo-Hookean model to describe shear deformations in elastomers [7,8]. The Fung model was developed initially to capture the response of tissues with a high content of collagen fibres, such as skin and arterial walls [9,10]. This model exhibits a typical dramatic strain-stiffening response in uniaxial loading characterizing the extension of stiff crimpled collagen fibres. The Gent model further penalizes this extension by limiting the strain to a finite value, similar to worm-like chain models used in polymer physics [11]. As such, both the Gent and Fung models are suitable for tissues that derive their elasticity from a mixture of a soft elastic matrix and stiff fibres. However, aggregates of cells found in brain and in fat tissues are approximately equiaxed in structure with a large lipid content, and this accounts for their almost isotropic, incompressible properties, which likely originate from their cellular structure rather than fibres [4,12–15] (figure 1 and appendix A).

In this study, we provide a set of model examples for which we compare the values of the shear modulus under increasing compression or tension with experimental data for brain and fat tissues. The viscoelasticity of brain and adipose tissues was measured following the protocol described in [4]. The dynamic shear storage modulus *G*′ was measured as a function of time for increasing tensile or compressive strain (from 0% to 40%). Details are given in appendix A.

A hyperelastic constitutive material has a unique stress–strain relationship, independent of strain rate. However, the stress–strain response for viscoelastic materials changes with strain rate, and a strain–energy density function does not exist for these materials. Nonetheless, for many soft tissues, the shape of the nonlinear stress–strain curve is typically invariant with respect to strain rate. In this case, at fixed strain rate, the shear modulus may be captured by a nonlinear hyperelastic model (an example of this approach for fat tissues can be found in [16]).

The usual practice in constitutive modelling is to fit uniaxial data obtained under controlled compression, and less commonly under tension, with standard material models. This is due to the general experimental (and partly analytical) limitations to carry out proper assessment of stresses and deformations in multiple loading situations. This approach has been particularly successful for tissues that operate mostly under axial loading conditions, such as tendons or ligaments. However, soft tissues such as fat or brain, operate in highly varying and complex loading environments and exhibit responses that cannot be easily modelled by such an approach. The elasticity of these materials can be probed by subjecting samples to multiple loadings and, indeed, our study of the shear modulus under combined stretch and shear demonstrates that different constitutive models behave very differently under combined deformation, even though some of them may respond very similarly in axial deformation alone.

From our numerical results, we infer that the shear modulus for the Mooney material is too small compared with the experimental values at similar strains, but appropriate Ogden models are found which are in excellent agreement with the experiments, and thus conditions (i)–(ii) are satisfied by the corresponding shear modulus. The newly identified models are robust and suitable for use in large-scale finite-element computations. Furthermore, for the Mooney, Fung, Gent and Ogden models analysed here, the elastic modulus increases or remains almost constant as compression increases, and therefore, condition (iii) is satisfied numerically, whereas for the neo-Hookean material, this modulus decreases under increasing compression.

For the hyperelastic models under consideration, the associated strain energy functions and their ability to satisfy the conditions (i)–(iii), either theoretically or numerically, are summarized in table 1. The numerical results are shown at a glance in figures 3 and 4.

## 2. Nonlinear elastic modulus and shear modulus relations

The homogeneous (affine) deformations analysed here are universal and controllable in the sense that they can be maintained in every homogeneous, incompressible, isotropic, elastic material by application of suitable surface tractions [22–26]. If the material is described by a strain energy function the associated Cauchy (true) stress has the Rivlin–Ericksen representation
where *p* is the arbitrary hydrostatic pressure, **B** is the left Cauchy–Green strain tensor with the principal invariants *I*_{1}, *I*_{2}, *I*_{3} and
are the material response coefficients. Equivalently, in terms of the principal stretches *λ*_{1}, *λ*_{2}, *λ*_{3}

Henceforth, we assume that these material responses are consistent with the Baker–Ericksen (BE) inequalities stating that *the greater principal stress occurs in the direction of the greater principal stretch*, and the pressure–compression (PC) inequalities stating that *each principal stress is a pressure or a tension according as the corresponding principal stretch is a contraction or an elongation* [27–29].

### 2.1. The elastic modulus in finite tension or compression

We first consider a unit cube of incompressible hyperelastic material subject to the uniaxial tension or compression in the second direction
2.1where (*x*, *y*, *z*) and (*X*, *Y*, *Z*) are the Cartesian coordinates for the current and the reference configuration, respectively, and *a* > 1 (tension) or 0 < *a* < 1 (compression) is constant.

For the deformation (2.1), the left Cauchy–Green strain tensor takes the form
and the non-zero components of the associated Cauchy stress are
We define the *nonlinear elastic modulus* in the second direction as the ratio between the Cauchy (true) stress *σ _{yy}* and the logarithmic (true) strain (the sum of all the small strain increments) [24, p. 118]
2.2

If *σ _{xx}* =

*σ*= 0, then (2.2) takes the form 2.3and if

_{zz}*σ*=

_{xx}*σ*≠ 0, then by the PC inequalities,

_{zz}*σ*< 0 when 1/

_{zz}*a*< 1, and

*σ*> 0 when 1/

_{zz}*a*> 1, hence 2.4

### 2.2. The shear modulus for finite shear superposed on axial stretch

We further examine a unit cube material sample deformed by the combined stretch and shear
2.5where (*x*, *y*, *z*) and (*X*, *Y*, *Z*) are the Cartesian coordinates for the deformed and the reference configuration, respectively, and *a* and *k* are positive constants representing the axial stretch and the shear parameter, respectively (figure 2).

For the deformation (2.5), the left Cauchy–Green strain tensor takes the form and the non-zero components of the associated Cauchy stress are For the deformed cube, the shear strain on the inclined faces and the associated shear traction are, respectively

We define the *nonlinear shear modulus* as the ratio between the shear traction *σ _{t}* and the logarithmic shear strain ln(

*B*+ 1), i.e. 2.6

_{t}Then, the shear modulus (2.6) is independent of the hydrostatic pressure –*p* and is positive if and only if *β*_{1} − *β*_{−1}/*a* > 0.

When the shear strain is small, the shear modulus (2.6) takes the form 2.7where

Assuming that *σ _{zz}* = 0, we also define
and obtain
2.8Therefore, as the axial stretch

*a*increases, the magnitude of the normal force

*N*

_{0}relative to the shear modulus

*μ*

_{0}also increases. This is a

*universal relation*[30], which holds independently of the material responses

*β*

_{1}and

*β*

_{−1}, and is analogous to

*Rivlin's formula*for a cylinder deformed by combined stretch torsion [26, p. 192]. Then, by (2.3) and (2.8) 2.9i.e. the ratio between the elastic modulus

*E*and the shear modulus

*μ*

_{0}is also independent of the material parameters, and

*E*/

*μ*

_{0}→ 3 as

*a*→ 1.

If *σ*_{zz} ≠ 0, then by the PC inequalities, *σ _{zz}* < 0 when 1/

*a*< 1, and

*σ*> 0 when 1/

_{zz}*a*> 1. In this case hence Then, by (2.4)

### 2.3. The behaviour of nonlinear hyperelastic models

For the hyperelastic materials listed in table 1, we examine the elastic modulus (2.3) and the shear modulus (2.7) as the magnitude of the compressive or tensile strain *b* = ln *a* increases. In view of the subsequent comparison with experimental data, we restrict our attention to the case when

— For the neo-Hookean model, the shear modulus (2.7) is equal to and is independent of strain.

*Hence, condition*(*ii*)*is satisfied, but not*(*i*).Applying (2.9) the corresponding elastic modulus takes the form 2.10and increases as increases.

*Thus, condition*(*iii*)*is not satisfied*.— For the Mooney–Rivlin model, the shear modulus takes the form and, if

*C*_{1}> 0 and*C*_{2}> 0, then this modulus decreases as*b*increases.*Hence, conditions*(*i*)*and*(*ii*)*are both valid*.— By (2.9), the elastic modulus is equal to 2.11and, if then this modulus decreases as increases.

*Thus, condition*(*iii*)*is also valid*.— For the Fung model: and, if

*C*> 0 and*α*> 0, then this modulus decreases as increases and increases as increases.*Hence, condition*(*i*)*is satisfied, but not*(*ii*).For this model, by (2.9), the elastic modulus is 2.12and there exists , such that this modulus decreases as increases and increases as increases.

*Thus, condition*(*iii*)*is not satisfied*.— Similarly, for the Gent model: and, if

*C*> 0 and*β*> 0, then this modulus decreases as increases and increases as increases.*Hence, condition*(*i*)*is satisfied, but not*(*ii*).By (2.9), the corresponding elastic modulus is 2.13and there exists , such that this modulus decreases as increases and increases as increases.

*Hence, condition*(*iii*)*is not satisfied*.— For the Ogden model: where For this model, a general conclusion about the monotonicity of the shear modulus cannot be drawn, and particular cases need to be examined individually. We do this in §3 where hyperelastic models are treated numerically.

## 3. Numerical results

Here, we compare the mechanical performance of the neo-Hookean, Mooney, Fung, Gent and Ogden materials described above when fitted to available experimental data for the shear modulus of brain and fat tissues. According to the experimental measurements, the shear modulus increases strongly under increasing compression, whereas in tension, it remains almost constant or decreases slightly at first, then begins to increase, but much less than in compression. In particular, for brain tissue, the shear modulus is essentially constant up to 10% tensile strain, whereas for lean and obese fat, this modulus appears almost constant up to 30% and 20% tensile strain, respectively. The experimental data for brain and fat tissues are marked by the (red) circles in the plots shown in figures 3*a* and 4*a*, respectively.

The values of the constant parameters for the hyperelastic models fitted to brain and fat data are recorded in tables 2–3 and tables 4–5, respectively. For the neo-Hookean, Mooney–Rivlin, Fung and Gent models, all constant parameters were fitted. For the Ogden models, the non-zero coefficients *C _{p}* were fitted, whereas the corresponding exponents

*m*were fixed. The fitting of the material parameters was performed using a nonlinear least-squares procedure implemented in Matlab (

_{p}*lsqnonlin.m*). By this procedure, the following (unconstrained) minimization problem was solved where

**c**= (

*c*

_{1},

*c*

_{2}, …

*c*) are the constant material parameters to be identified, (

_{m}*b*,

_{i}*μ*) are the pairs of data for the compressive or tensile strain and the shear modulus, respectively, and

_{i}*G*(

*k*,

_{i}*b*;

_{i}**c**) =

*μ*(

*k*,

_{i}*a*) is the shear modulus defined by (2.6), such that is the stretch parameter, and the shear strain is constant and small, viz. 0.02 for brain and 0.035 for fat tissues, so we set

_{i}*k*= 0.02

_{i}*a*and

_{i}*k*= 0.035

_{i}*a*, respectively, for all

_{i}*i*= 1, …,

*n*.

In order to assess the accuracy with which the models capture the mechanical behaviour measured by the experiments, for each model, the relative error of the shear moduli to the given data was also computed, as follows 3.1

For the neo-Hookean, Mooney, Fung and Gent models with constant parameters as indicated in table 2, the shear moduli at 2% shear combined with up to 40% compression or tension, and their relative errors (3.1) are plotted in figure 3*a*,*b*, respectively. Because the shear strain is small, the shear modulus *μ*_{0} defined by (2.7) is capable of predicting theoretically the corresponding mechanical behaviour of these models under the combined deformation. For these models, we further compute the elastic modulus *E* defined by (2.3), and plot its values normalized to those at 5% compression in figure 3*c*. Numerically

— For the neo-Hookean material, the computed shear modulus

*μ*defined by (2.6) is virtually constant,*hence condition*(*ii*)*is valid, but*(*i*)*is not*. For this material also, the relative values of the elastic modulus (2.10) plotted in figure 3*c*decrease when compression increases,*hence condition*(*iii*)*is not valid*. These results are all in agreement with the theoretical findings for the neo-Hookean model.— For the Mooney–Rivlin material, the shear modulus

*μ*increases as compression increases and decreases as tension increases,*thus conditions*(*i*)*and*(*ii*)*are both satisfied*. From the relative values of the elastic modulus plotted in figure 3*c*, we also see that, for this material, the elastic modulus (2.11) increases under increasing compression,*i.e. condition*(*iii*)*is also satisfied*. These results are again in agreement with the theoretical findings for the Mooney model. Unfortunately, the numerical values of the shear modulus attained by this model are much smaller than those required by the experimental results for brain tissue, as shown by the large relative error estimates, hence Mooney materials, which have proved excellent in describing elastomers and other materials with entropic elasticity, are inadequate for the modelling of this tissue.— For the Fung and the Gent materials, the respective shear moduli

*μ*increase as compression increases, i.e.*condition*(*i*)*is satisfied*, but because they also increase almost as fast in tension as in compression,*condition*(*ii*)*is not satisfied*. Moreover, the corresponding relative errors increase rapidly as either compression or tension increases, hence, these materials do not capture the required physical behaviour in either of these deformations. For these models also, the monotonicity of the associated elastic modulus (2.12) and (2.13) changes, albeit slowly, so that the computed modulus remains almost constant before it increases as compression increases.*Hence*,*condition*(*iii*)*is, in fact, satisfied numerically*.

As the neo-Hookean, Mooney, Fung and Gent models fail to agree with the experimental results for brain tissue under combined stretch and shear, and similar results are shown to hold experimentally for adipose tissue, we illustrate numerically the behaviour of these material models in rapport to the brain data, but take these models no farther when modelling fat tissues (table 6).

— For brain and fat tissues, we further determine six different Ogden-type models, with the associated constant parameters recorded in tables 3–5. In tables 3–5, the Ogden

models have_{N}*N*non-zero coefficients*C*, whereas the associated exponents_{p}*m*are fixed. For these models,_{p}*conditions*(*i*)*and*(*ii*)*are both valid*. See also figures 3*a*and 4*a*. The relative errors recorded in tables 7–9 further suggest that Ogden_{7}and Ogden_{8}are the most successful in approximating the experimental data. Obviously, these last models contain a large number of parameters and are likely to over-fit the data. The purpose of including these models is to demonstrate that such a family of models is adequate to capture the mechanical responses of the biological tissues under investigation. See also figures 3*b*and 4*b*. From the associated relative elastic modulus plotted in figures 3*c*and 4*c*, we also see that this modulus increases under increasing compression,*hence condition*(*iii*)*is also valid*. For all models, smaller relative values for the elastic modulus may be obtained when this modulus is defined by (2.4).

As explained above, for brain and fat tissues, it was observed experimentally that the shear modulus increases sharply under increasing compressive strain, but not under tensile strain, whereas the apparent elastic modulus increases or remains almost constant when compressive strain increases. The macroscopic, centimetre-scale, samples are heterogeneous on a smaller length scale, because they are a mix of grey and white matter with boundaries between, but this heterogeneity does not dominate the rheological response, because grey and white matter do not differ strongly in stiffness, and the macroscopic viscoelastic response does not depend on precisely how the sample is cut or how large it is. In addition, a characteristic of these tissues is that they exhibit a predominantly isotropic incompressible behaviour. Here, we compare the behaviours of several nonlinear hyperelastic models, both theoretically and numerically, and test our results against available experimental data. Our analysis shows that neo-Hookean, Mooney–Rivlin, Fung and Gent models, which have been successfully employed to date in the modelling of rubber and of other man-made or natural materials, are inadequate to model brain and fat tissues. Instead, for these tissues, Ogden models, with four, six and eight coefficients, respectively, are found which are in excellent agreement with the experiments. The newly identified models can be easily implemented in finite-element codes.

## 4. Conclusion

Biological tissues offer a great diversity of mechanical responses when subject to loads. Often, such behaviours appear counterintuitive as our intuition has been forged by centuries of studies of engineering material often treated in the limit of small strains. It is then tempting to conclude that classical continuum mechanics is not suitable for modelling biological materials. However, aside from the very few classical models used indiscriminately for both rubbers and biological tissues, there is a vast pool of potential models that have yet to be explored, understood and classified.

Our approach when confronted with a new constitutive phenomenon consists of two steps: first, based on experimental evidence, classify qualitative responses that a model ought to reproduce by an analytical study of relevant deformations. Second, for the models that pass the first sift, find suitable candidates in quantitative agreement with the data. As presented here, this approach was successful for the analysis of the response of brain and fat tissues, generating in the process models that are capable of predicting some of the key elastic properties underpinning the extraordinary mechanical performance of these tissues, and which can be integrated into large-scale computational framework. Far for claiming that the models presented here are universal models for brain and fat tissues, we demonstrate that a systematic approach in the framework of nonlinear elasticity based on experimental data provides phenomenological models that can be used to explore the large-scale response of tissues and organs. We have listed several models that fit the data increasingly well at the expense of an increase in the number of parameters. We leave it to the practitioners to decide, based on the problem at hand and the range of deformation being studied, which model to use.

Our enquiry also suggests that the microscopic processes that generate macroscopic elastic responses in these tissues are different to the ones found in other soft tissues or elastomers. Clearly, there is a need for better understanding of the mechanics of very soft tissues, particularly of the brain, which is currently under intense study by researchers in both biophysics and computational mechanics, and adipose tissue, which is a growing area of investigation in clinical research.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material.

## Authors' contributions

A.G. and P.A.J. conceived and designed the study; L.A.M. and A.G. carried out the hyperelastic analysis and drafted the manuscript; L.K.C. and P.A.J. carried out the experimental work and helped draft the manuscript. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

A.G. is a Wolfson/Royal Society Merit Award Holder and acknowledges support by the Marie Curie European Reintegration grant no. BKRVRG0. The support for L.A.M. by the Engineering and Physical Sciences Research Council of Great Britain under research grant no. EP/M011992/1 is gratefully acknowledged.

## Acknowledgements

We thank the hospitality of the NSF-Mathematical Biosciences Institute in Ohio for the workshop where the ideas presented here were first discussed, as well as key discussions with Prof. Ellen Kuhl of Stanford University.

## Appendix A. Experimental method

Brain and lean adipose tissues were harvested from adult male wild-type C57BL/6 mice and obese adipose tissue was collected from genetically obese ob/ob mice (The Jackson Laboratory, Bar Harbor, ME). Fresh tissues were stored in Dulbecco's modified Eagle's medium (DMEM, Gibco, Grand Island, NY) and tested within a maximum of 3 h after sacrifice. Macroscopic rheometry using a Rheometrics fluids spectrometer III strain-controlled rheometer (Rheometrics, Piscataway, NY) fitted with 8 mm diameter parallel plates was used to measure the viscoelasticity of tissue samples using methods previously described in [4]. During testing, tissues were kept hydrated by surrounding the sample with phosphate-buffered saline (PBS). Control experiments showed that the use of DMEM and PBS increased tissue weight only slightly—on average 4% for brain, 2% for fat.

Test samples were cut into disc-shaped samples using an 8 mm diameter stainless steel punch, maintaining the cerebral hemispheres, but removing the cerebellum, olfactory bulb, pons and medulla (figure 5). Note that an 8 mm diameter disc in the brain will necessarily contain a mixture of different tissues and structures. However, structures such as ventricles and vessels constitute a small percentage of the volume (2%) and are not arranged in a highly oriented manner. The bulk of the tissues consists of cells packed together and surrounded by soft matrices, which in the brain are flexible polysaccharides and so very likely to be close to isotropic (as found repeatedly in other experiments [1]). Because the length scale of the individual tissue components is much smaller than the length scale of the bulk tissue, the protocol provides an average mechanical response suitable for comparison with continuum models. If swelling occurs and we assume that it deforms the sample equally in all directions, a typical 4% increase in volume yields 1% increase in height. A 1% change in height and the effect on *G*′ are relatively small compared with the deformations the samples are subjected to during measurement.

To avoid sample slippage during shear deformation and to perform uniaxial extension, fibrin gel was used to glue the sample to the rheometer plate. Fibrin gel was prepared by mixing equal volumes of 28 mg ml^{−1} salmon fibrinogen solution with 125 U ml^{−1} of thrombin (Sea Run Holdings, Freeport, ME) directly on the lower rheometer plate, and the sample was immediately positioned. Subsequently, a thin layer of fibrin gel was pipetted onto the upper surface of the tissue, and the top plate was lowered until a positive normal force (1 g) was measured by a force transducer. Control experiments showed that addition of fibrin glue did not affect the viscoelastic properties of the tissue samples.

The dynamic shear storage modulus *G*′ was measured as a function of time (brain: 2% oscillatory shear strain, 2 rad s^{−1} frequency; fat: 3.5% shear strain, 2.5 rad s^{−1} frequency) and increasing tensile/compressive strain (0–40%; figure 5). For both brain [4] and lean fat, recoverable deformation was observed when the compressed sample was returned to its initial height and the shear modulus *G*′ returned to its uncompressed value (figure 6). Obese fat largely recovers, but residual stresses appear to remain after 40% compression. This may be indicative of tissue damage at higher levels of compression, but does not negate the compression-stiffening phenomenon observed. During each incremental compression, *G*′ demonstrated a relaxation response. Therefore, for the purposes of model fitting, *G*′ values after 100 s of relaxation were used. In tension, the following correction was applied to account for the decrease in cross-sectional area during testing under the assumption that volume is conserved, where *G*′ is the storage modulus and *λ* is the axial strain. To determine the linear viscoelastic region for each tissue, the shear moduli *G*′ of brain and lean fat were measured with respect to increasing shear strain up to 10% (*n* = 3 each, figure 7). The experimentally measurable linear viscoelastic region for brain and fat were determined to be 0.15–2.5% and 0.15–4% strain, respectively.

- Received May 30, 2015.
- Accepted August 20, 2015.

- © 2015 The Authors.

Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.