In this paper, we review critically some basic models derived from the statistics of long-chain molecules and then discuss the status of such models within the three-dimensional nonlinear theory of elasticity. We draw attention to some deficiencies of certain worm-like chain (WLC) models when viewed within the three-dimensional continuum mechanics framework. Modifications of the WLC models motivated by such considerations and avoiding these deficiencies are then discussed and shown to correspond well with data generated by the exact WLC model.
Rubber elasticity has one of the longest and distinguished histories in all of polymer science (see, for example, the discussion in Saccomandi & Ogden 2004). However, the complex molecular nature of rubbers (and other elastomeric materials) and the need for the use of nonlinear theories to describe their mechanical behaviour have made difficult the development of realistic mathematical models. The theory of incompressible, isotropic hyperelasticity is a starting point for the modelling of elastomeric materials, and within this framework two main directions can be followed in developing specific constitutive models. These are the phenomenological approach based on continuum mechanics arguments (Saccomandi 2004), and the statistical approach based on molecular concepts (Erman 2004).
In following the molecular approach, the first step is to derive an elasticity model for a single polymeric chain within the rubber network and then to introduce a suitable averaging procedure to obtain a macroscopic constitutive equation for the material response. For example, in the elementary molecular theory of networks it is postulated that the elastic free energy of a network is equal to the sum of the elastic free energies of the individual chains (Erman 2004). On the other hand, the main tools of the phenomenological approach are the axioms of continuum mechanics together with a systematic and rigorous use of the restrictions imposed by thermodynamics, the use of representation formulae based on symmetry arguments and a priori restrictions induced by various mathematical requirements (Saccomandi 2004).
Rubber-like elasticity is also fundamental in the study of biological proteins and materials that are characterized by high resilience, large strains and low stiffness (Gosline et al. 2002). This is, for example, the case in the pioneering paper of Fung (1967) on mesentery tissue. In his paper, which may be considered as a starting point for the modern theory of cardiovascular solid mechanics (Humphrey 2002), an empirical relation for the highly nonlinear uniaxial stress-stretch behaviour of the tissue is proposed. This relation is the origin of the celebrated exponential Fung model for soft tissue mechanics (e.g. Humphrey 2002). However, this one-dimensional empirical relation is not consistent with the equations of nonlinear elasticity since the associated uniaxial stress does not contain the factor λ−1 that arises for simple tension in the general theory (λ being the stretch in the direction of extension). Thus, without modification, it cannot be extended to the general three-dimensional setting that is required to ensure its consistency within the more general framework. What is today referred to as the Fung model in this general setting was introduced, in the case of isotropy, by Demiray (1972). In recent years, there has been a steady stream of theoretical and experimental research aimed at constructing models compatible with the laws of continuum mechanics but also related, where possible, to the molecular approach (e.g. Holzapfel & Ogden 2003).
The situation that we have described above for soft tissue mechanics is similar to that arising in consideration of single-molecule studies within the context of DNA mechanics and the mechanics of other biological macromolecules and proteins. Stretching of individual biomolecules is now achieved by a variety of techniques that allow measurement of forces from 10 femtonewtons to hundreds of piconewtons (Bao & Suresh 2003). The entropic elasticity of these biomaterials is usually described by using one-dimensional models based on statistical mechanics expressions of the entropy that idealize each macromolecule as a freely jointed chain (FJC) or a worm-like chain (WLC).1 This is the case for DNA (Strick et al. 2003) and also for some structural proteins such as elastin (Urry et al. 2002) and collagen (Sun et al. 2004). A detailed introduction to the WLC may be found in appendix G of the book by Flory (1988).
The idea of using a general three-dimensional model based on a chain network in which the basic building block is the WLC has recently been proposed by Bischoff et al. (2002a,b) and Kuhl et al. (2004). The aim of the present paper is to build upon this work by considering a three-dimensional version of the WLC model within the context of the continuum theory of nonlinear elasticity.
In §2, we begin by reviewing critically some basic models derived from the statistics of long-chain molecules with particular attention to problems that may arise in fitting the models to the data. This is followed by discussion of the status of such models within the three-dimensional nonlinear theory of elasticity. In §3, we draw attention to deficiencies, from the continuum mechanics perspective, of a specific WLC model, with particular reference to simple tension and simple shear. Sections 4 and 5 are concerned with modifications of the Gent model of rubber elasticity (Gent 1996) that are motivated by the above considerations and avoid the deficiencies encountered in the WLC model. Some closing remarks are contained in §6.
2. Basic equations
(a) Some long-chain molecular models
An introduction to the statistical properties of long, flexible objects such as polymer chains may be found in, for example, the book by de Gennes (1979). In the freely jointed chain model, a molecule is considered to be made up of rigid orientationally independent chains. Alignment of segments by a tension, denoted f, is described by using non-Gaussian statistics with the relative stretch r/L as a variable, where r is the end-to-end distance of the chain and L its contour length. Specifically,(2.1)where k is the Boltzmann constant, T the absolute temperature, N the numbers of freely jointed segments (each of length l, so that L=Nl), and 0<r≤L. The symbol −1 stands for the inverse of the Langevin function, which is defined by(2.2)This is a classical model from which various three-dimensional micromechanical network chain models have been derived (Boyce & Arruda 2000).
In the WLC model, the chain is treated as a flexible beam that curves as a result of thermal fluctuations. An excellent treatment of the WLC may be found in the paper by Bouchiat et al. (1999), where this model is derived by using modern methods of statistical physics. Because, in the case of the WLC, an analytically exact expression for the force in a single chain is not available, it is usual to consider an interpolation formula of the form(2.3)as introduced by Marko & Siggia (1995), where A, with l≤A≤L, is a parameter referred to as the persistence length of the chain. A more complex form than (2.3) is also given in Bouchiat et al. (1999), where an amended version is proposed that takes account of some enthalpic corrections for large forces. This is written(2.4)where αi, i=2, …, 7, are phenomenological parameters that account for the enthalpic corrections.
Here, we use the term ‘WLC model’ for these interpolation formulae instead of for the exact model, and it is therefore important to point out that our comments are relative to the interpolation formulae rather than the real WLC model.
At this point some comments on the formulae (2.1) and (2.3) are called for, for which purpose we introduce the extension ratio z=r/L, such that 0<z≤1. Equation (2.3) may be rewritten as(2.5)This is the ratio of a third- and a second-order polynomial and has a singularity at z=1. Also f(0)=0. In (2.4), a polynomial term is added to the basic function (2.3) in order to improve the overall fitting of the model to the data, and we now discuss this ‘improvement’ in detail. Our aim is to show the importance of the singularity in (2.3) and for this reason we introduce the reduced WLC function(2.6)i.e. equation (2.5) with the term 4z omitted.
In figure 1, log plots of both (2.5) and (2.6) are shown for comparison with the ‘data’ taken from Bouchiat et al. (1999). These are data generated numerically from the exact WLC model. The least squares fitting problem can be expressed in general as minp F(p), where the squared two-norm of residues (the ‘objective function’) is defined by(xi, yi), i=1, …, m, are the data points and p=[p1, …, pn] represents the vector of parameters to be identified, f(p,xi) is the force predicted by the constitutive model for the value xi and yi is the corresponding experimental value of the force. If p* is the optimum set of parameters obtained by the above fitting procedure the corresponding residual is S=F(p*) and the percentage relative error is defined byIn order to evaluate the ‘goodness of fit’ it is also usual to calculate the root mean square error (RMSE), here defined aswhere m−n is the number of degrees of freedom of the optimization problem and fref is the mean of the yi data values. Note that, in contrast to S, ϵ is dimensionless. In the following, we also report this value for each example.
When we deal with equations (2.5) and (2.6), the only parameter to be identified is kT/4A and in this case the problem can easily be solved as a linear problem. To deal with the nonlinear least squares problems in the following sections we use the lsqcurvefit routine in Matlab 6.5 by selecting the Trust Region Method.
In the fitting reported in figure 1, we record that for (2.5) we have kT/4A=0.2486 and that the residual S is 2.3295 and the RMS error is ϵ=0.0106, while for (2.6) the corresponding values are kT/4A=0.2506, S=1.2018 and ϵ=0.0076. Thus, the residual S and the RMSE ϵ are less for the reduced form (2.6), although this is not apparent from inspection of the curves since the fit is clearly better for (2.5) in the lower range of extensions. The differences in S and ϵ are accounted for by the much better fit of (2.6) at larger extensions. This is clear from the numerical values but is not apparent from the figure because of the large gradients in this region. Our analysis shows that the function in (2.6) needs to be improved but that ‘improvement’ by addition of the term 4z, introduced to fit the data for small extensions, is not good. Indeed, the quality of the fitting deteriorates at the higher values of the extension although a better agreement for small extensions is achieved. For this reason, the formula (2.4) was proposed by Bouchiat et al. (1999). We note, however, that the large number of constitutive parameters introduced in (2.4) may be problematic. Indeed, in Ogden et al. (2004) it was shown that increasing the number of polynomial terms in phenomenological models may introduce numerical instabilities and non-uniqueness of the optimal sets of parameters. The problem of multiple solutions has not been apparent previously in the present context because DNA elongation data have been based on the formula (2.4) with the coefficients αi set to the values identified by Bouchiat et al. (1999). In fact, (2.4) is normally used in a nonlinear fitting procedure to identify only the contour length L and the persistence length A of actual data. We deal with this point in a forthcoming paper (Ogden et al. submitted). It should be emphasized that it is the term with the singularity that is responsible for the very good representation of the non-Gaussian behaviour of the data at larger extensions.
We end this section by recalling that the elastic energy function W associated with (2.5) is(2.7)where W0 denotes a constant energy with respect to a reference level. This is a one-dimensional strain energy. In §2b, we consider the possibility of extending this model to three dimensions.
(b) Nonlinear elasticity and molecular models
It is important to point out that the expressions for the tensile force in equations (2.1), (2.3) and (2.4) are not compatible with the three-dimensional theory of elasticity. This is not just an academic consideration since there are many situations where single molecule experiments are performed involving torsional deformations (e.g. Strick et al. 2003), and torsion requires three-dimensional analysis. Within the framework of the phenomenological theory of nonlinear elasticity the constitutive equations for the strain-energy function of an isotropic material are based on the choice of a scalar function of a suitable set of deformation invariants or, equivalently, on the principal stretches of the deformation. The usual invariants adopted are the principal invariants of the right Cauchy–Green deformation tensor C=FTF, where F is the deformation gradient. These are defined by(2.8)In terms of the principal stretches λ1, λ2, λ3, we have(2.9)When the internal constraint of incompressibility is in force all admissible deformations must satisfy the restriction . A basic result from the general theory of nonlinear elasticity (e.g. Ogden 1984) is that for an incompressible isotropic elastic material characterized by a strain energy W(I1, I2), the Cauchy stress tensor T has the form(2.10)where I is the identity tensor, B=FFT is the left Cauchy–Green deformation tensor and Wi=∂W/∂Ii, i=1, 2.
It is clear that any possible incompressible isotropic strain energy derived by molecular arguments must be expressed as a function of the first two invariants in (2.9), and the corresponding Cauchy stress tensor must be compatible with formulae such as (2.10). An important example derived by using molecular arguments starting from the non-Gaussian FJC theory is that due to Arruda & Boyce (1993), while Beatty (2003) has provided continuum formulations of the various strain energies that are based on the FJC model.
Consider next a single perfectly flexible and freely jointed chain whose undeformed mean length is assumed to be . The relative chain stretch λr is the ratio of the current chain vector length r to its fully extended limiting length Nl. Thus,(2.11)where λchain=r/r0. Note that λr∈[N−1/2,1]. Now suppose that the chain vector r0=(X, Y, Z) in the reference configuration is directed along the line X=Y=Z through the origin of a rectangular Cartesian system. This has length and its direction is . If the network is subjected to an affine deformation in which the rectangular Cartesian axes coincide with the principal axes then the squared stretch of a chain initially oriented in an arbitrary referential direction m is given by(2.12)which specializes to I1/3 for the specific m above. The same result holds for a chain whose end point is oriented initially according to −X=Y=Z or −X=−Y=Z or X=−Y=Z. For the Arruda–Boyce eight-chain model, we then have(2.13)
Beatty (2003) showed that some of the constitutive equations proposed on the basis of the FJC model lead to a Cauchy stress tensor in the simple invariant form(2.14)where(2.15)For a single FJC, we have(2.16)
To avoid the use of the inverse Langevin function, it is usual to replace −1(λr) by a polynomial expansion (e.g. Boyce 1996), for which the crucial characteristic of the singularity in (2.16) is then lost. Another possibility is to replace the inverse Langevin function by its Padé approximation, as, for example, in Miehe et al. (2004). However, Horgan & Saccomandi (2002a) have shown that a very good global approximation is obtained by using the phenomenological model first proposed by Gent (1996). The Gent model has strain energy(2.17)where μ is the shear modulus for infinitesimal deformations and Jm(>0) is the limiting value of I1−3 associated with limiting chain extensibility. (The upper bound on I1 is consistent with that associated with (2.13) if we take Jm=3N–3.) In the limit, as the chain extensibility parameter tends to infinity (Jm→∞), (2.17) reduces to the classical neo-Hookean form(2.18)
The model (2.17) is a fundamental one because it is associated with the simplest rational approximation of the response function in the Cauchy stress representation formula for an incompressible isotropic elastic material for which there is a fixed singularity (Horgan & Saccomandi 2003). The Cauchy stress associated with (2.17) is given by(2.19)so that the stress has a singularity as I1→Jm+3. Note that this is a special case of (2.14).
3. An eight-chain WLC-based model
Kuhl et al. (2004), following Bischoff (2002a,b), have developed a three-dimensional version of the WLC model by considering a three-dimensional strain energy associated with (2.3). This may be obtained by establishing first a connection between r and the invariants in (2.9) and then relations similar to (2.11) and (2.13), as for an eight-chain WLC model. In remark 1 of Kuhl et al. (2004), the relationis considered, where a is a constant parameter that may be chosen so that r/L coincides with (2.13). It should be noted that the range of values of r/L in the one-dimensional WLC models is (0,1] for the extension ratio, whereas the stretch in the direction of extension is greater than one.
From equation (2.7), we obtain(3.1)Now, for a fixed temperature, the persistence length A is a constant proportional to the bending stiffness of a chain segment. Then, when N→∞ in (3.1), we find that this strain energy blows-up, which is a pathological and un-physical result. It might be expected that for N→∞ we recover the neo-Hookean energy function or at least a ‘classical’ strain energy that is well defined in this asymptotic limit (Horgan & Saccomandi 2003). It does not seem possible to circumvent this problem for this model, even if A is treated as a phenomenological parameter that depends on N.
In what follows it is convenient to introduce the notation J0=3N and, for curve fitting purposes, to treat kTNl/4A as a single parameter. Then, we re-write (3.1) as(3.2)from which we deduce that(3.3)which will be required in §3a.
In the paper by Kuhl et al. (2004), the basic model (3.1) is considered in a more general framework in which anisotropy, compressibility and growth are included, and their paper focuses more on the description of remodelling of biological tissue than on a detailed study of the mechanical properties of the strain energy associated with the WLC model. In the subsequent paper by Kuhl et al. (in press), the model derived within their framework is analysed mathematically with a view to establishing various convexity properties.
Here, we study some basic mechanical properties of (3.1) and, to this end, we start by considering the mechanical response of this model for some simple homogeneous deformations.
(a) Simple tension
For an incompressible material, the principal stretches in simple tension are (e.g. Ogden 1984)(3.4)where λ is the stretch in the direction of extension. It follows that(3.5)and the principal components of the Cauchy stress tensor (2.10) are given by(3.6)The requirement that the lateral surfaces of the specimen undergoing the simple extension are traction free, t2=t3=0, yields(3.7)and the tensile force f per unit of undeformed cross-sectional area necessary to achieve the stretch is given by(3.8)We emphasize that here and subsequently f is, therefore, the nominal (or engineering) stress, whereas in the one-dimensional models it is the tensile force.
It is interesting to compare (3.9) with the tensile stress associated with the Gent material (2.17), namely(3.11)In the case of (3.11), the limit Jm→∞ enables us to recover the tensile force(3.12)associated with the Gaussian (neo-Hookean) model (2.18).
Again it is possible to propose a reduced form of the WLC model by omitting the first and third terms in (3.9), which then reduces to(3.13)
Figure 2 shows a comparison of (3.9) and (3.13), which parallels that shown in figure 1 in respect of (2.5) and (2.6). In the case of (3.9), we have kTNl/4A=0.0875 and J0=5 with a residual S=0.041 and RMSE ϵ=0.0014, while for (3.13) kTNl/4A=0.0878, J0=5 and the residual is S=0.1316 and ϵ=0.0026 (figure 2c shows the relative errors). There is, therefore, no essential difference between the reduced and the full versions of the three-dimensional WLC, which can be appreciated from the log plot in figure 2b. In contrast to figure 1, however, the full version has the smaller residual and RMSE. A better agreement between the theoretical model and the data is obtained than for the one-dimensional model. This is related, in particular, to the fact that since the one-dimensional tension functions proposed in order to give an explicit representation of the WLC were not derived in a rigorous way from the three-dimensional theory they do not contain the important factor (λ−λ−2).
In figure 3, a comparison between the one-dimensional formula (2.5), the three-dimensional WLC (3.9) and the Gent formula in (3.11) is shown, with the previous data shifted from the range (0,1) to the range (1,2), corresponding to the stretch λ. Note that for J0=5, the asymptotic value of the stretch calculated from (3.10) is λ*=2. These plots enable the differences between the models to be appreciated. We note, in particular, that the Gent model is unable to fit the data and that the stress derived from the three-dimensional model gives the better results for the considered data. For (2.5) and (3.9), the same material constants were used as in figures 1 and 2, respectively, while for (3.11) we found the values μ=1.2872, Jm=1.9188, with residual S=580.83 and RMSE ϵ=0.1709.
(b) Simple shear
Although simple shear deformations are not, per se, important in the context of biomolecular experiments, torsional deformations are important, and torsion is locally a simple shear. It is, therefore, of interest to investigate the response of the various models for a deformation that involves simple shear, which we do in this section. A simple shear deformation in the (X2, X3) plane is defined by(3.14)where κ>0 is the (constant) amount of shear (e.g. Ogden 1984). It follows that(3.15)and the squared principal stretches are computed as(3.16)
The shear stress is obtained from the specialization of (2.10) as(3.17)where the secant modulus of shear is defined as(3.18)The infinitesimal shear modulus is defined as(3.19)and we note that in the linear theory of incompressible isotropic elasticity this is the only constitutive parameter that arises.
In the case of (3.2), we have(3.20)and(3.21)
For the Gent model, we have(3.22)and therefore , a result that is independent of Jm∈[3,∞).
4. A WLC modification of the Gent model
Horgan & Saccomandi (2003) have shown that the model (2.17) proposed by Gent is indeed obtainable by a systematic procedure if rational approximations of the response functions for isotropic incompressible materials are considered. This result was obtained by generalizing the usual approximation method developed by Rivlin and others (e.g. Rivlin 1960) based on the expansion of the free energy as a truncated power series in the principal extensions ei, i=1, 2, 3, which are related to the principal stretches by λi=1+ei.
In terms of the principal extensions, the principal invariants (2.8) read(4.1)and the constraint is written as(4.2)In order to compute the power series it is convenient to introduce the equivalent set of invariants(4.3)from which it is easily shown that, when (4.2) holds,(4.4)Here, we have used the usual notation F=O(a) as a→0, if there exist constants δ>0 and ϵ>0 such that |F|<δ|a| whenever |a|<ϵ.
The method of Rivlin involves expanding the strain-energy function in the form(4.5)where Bαβ are suitable constitutive parameters and α and β may be integers or rational exponents. The technique proposed by Horgan & Saccomandi (2003) is based on approximating the response coefficients W1 and W2 in the representation formula (2.10) with the approximations based on rational functions instead of polynomials. The need for approximation is dictated by the fact that in experiments it is not the energy that is measured directly but quantities related to the first derivatives of the energy (forces or stress components). It is also, in general, necessary to consider supplementary conditions to ensure the existence of a strain-energy function, but such considerations are not needed here. We remark that it is well known that rational approximations may give more accurate results than polynomial approximations, in particular when singularities such as those associated with limiting chain extensibility models arise. Approximations using rational functions such as Padé approximants in this context have been described in Horgan & Saccomandi (2002b), where the relationships between the Gent model (2.17) and molecular models based on the inverse Langevin function are established. It is worth mentioning that in the recent review by Shaqfeh (2005) on the dynamics of the single-molecule DNA in flow, it was observed that the interpolation formula (2.3) is indeed a Padé approximant of the numerically simulated force obtained from the real WLC model.
(a) A 0/4 approximation
Our analysis of the interpolation form of the WLC suggests that a stronger form of singularity than that reported by Gent needs to be considered. We therefore take(4.6)where a1 and bi, i=1, 2, 3, are constants. This is a 0/4 approximation in the principal extensions. The constants in (4.6) are not arbitrary since they must satisfy some standard requirements of the nonlinear theory of elasticity. To reflect the finite chain length assumption a parameter J* is chosen so that as J1→J* both the stress and the strain energy associated with (4.6) become singular. It is also desirable that, as J*→∞, (4.6) and the associated strain energy should reduce to a standard polynomial form. Without loss of generality, the constants can be normalized by, for example, setting b1=1. We then have the approximation(4.7)where a and b are the two remaining (arbitrary) parameters.
The response coefficient (4.7) is associated with the strain-energy function(4.8)
In simple extension, the tensile stress calculated from (4.11) is(4.13)wherein J1=λ2+2λ−1−3. The factor J*–J1 has now been isolated to emphasize that the asymptote is fixed by the value of J* and is unaffected by the value of b. While μ has the interpretation as the infinitesimal shear modulus, the role of b is to control the rate at which the curves approach the asymptote and the flatness of the curve for small and moderate values of the stretch.
The fitting of this model to the data is considered in figure 4a as a log plot. Figure 4b shows the relative errors. We note that the infinitesimal shear modulus μ is very small (approximately 0.0052), while b=−0.4942 and J*=1.9803. The residual is S=0.1456 and the RMSE is ϵ=0.0028, so that it is clear that the fit for (4.13) is better overall than for the WLC model.
5. An improvement of the WLC modification
It is well known that generalized neo-Hookean strain-energy functions of the form W=W(I1) cannot describe completely the elastic properties of elastomeric materials. The criticisms of such models date back to Rivlin's investigation of natural rubber specimens (e.g. Saccomandi 2004), and have recently been reiterated by Horgan & Saccomandi (1999) and Wineman (2005) in considering the implications of a global universal relation for torsional deformations. Indeed, for generalized neo-Hookean materials the simplified representation formula for the Cauchy stress tensor involves only one response function, and this allows a new universal relation to be determined (i.e. a relation between the stress and strain components that is independent of the special form of the strain energy). This universal relation does not apply to a general incompressible isotropic energy function W=W(I1, I2). On comparison of experimental data with this special universal relation it can be deduced that the invariant I2 cannot be excluded if the energy function is to be regarded as realistic. For this reason, we examine here an enthalpic correction to the proposed strain energy involving I2.
For simplicity, we consider the two examples(5.1)and(5.2)where is defined in (4.11) and, in each case, C is an additional material constant. The first of these two models represents a generalization of the classical Mooney–Rivlin model, which is linear in I2, while the second is a generalization of the Gent–Thomas model (Gent & Thomas 1958). In Pucci & Saccomandi (2002), it was shown that for b=0 the model (5.2) allows the whole range of the experimental data in simple extension of a particular rubber to be described to an accuracy within 5%.
In figure 5, the stresses associated with (5.3) and (5.4) are fitted to the data on the same basis as for figure 4. Each model gives better results than (4.13). The values of the parameters are, for (5.3), μ=0.0020, b=−0.4977, J*=1.9912, C=0.0730, and, for (5.4), μ=0.0007, b=−0.4991, J*=1.9966, C=0.2880. The sums of residues are S=0.0179 and 0.0074, respectively; the RMS errors are ϵ=0.9998×10−3 and 0.6428×10−3, respectively.
6. Concluding remarks
In this paper we have considered the development of a phenomenological version of the WLC model. This is of interest for several reasons. First, because the WLC models used in the literature are in fact interpolation formulae of the real WLC model and therefore have a parallel phenomenological status to the continuum mechanics models. Second, the WLC models are essentially one-dimensional, whereas the models discussed in the present paper are three-dimensional and without the drawbacks that arise in the direct translation of the WLC classical interpolation formulae to three dimensions. Our methodology enables us to find three-dimensional phenomenological models that admit the accurate fitting of the ‘data’ generated numerically from the real WLC model. To appreciate how such models work for real experimental data, we consider finally an application to data from an actual single molecule experiment.
In figure 6, we show a set of data taken from Baumann et al. (2000) compared with the predictions of the model (4.13). These data have also been used recently by Rosa et al. (2005) to test a new interpolation formula for semi-flexible polymers. The data have been re-scaled by dividing the abscissa values by the first abscissa value so that the stretch λ can be used for the abscissa. Some general comments are appropriate at this point. Before approaching the vertical asymptote the data are very flat, which means that a very small value of the infinitesimal shear modulus is required, as was also the case with the ‘data’ used in §5. Moreover, from the data we can estimate the value of the limiting chain parameter. This information suggests possible starting values of the fitting parameters. In particular, we have linked J* to the location of the asymptote, while b is related to the flatness of the curve for f. This is an important issue since it is known that in nonlinear least squares fitting it is possible to have several optimal sets of parameters (Ogden et al. 2004) if random starting guesses in a wide range of values are considered. The parameters obtained here are μ=1.9466×10−5, b=−0.4995, J*=38.7247, the residual is S=0.8374 and the RMSE is ϵ=0.0587.
The approach that we have adopted allows the WLC model to be set within a three-dimensional continuum mechanical framework, and has been illustrated for some simple example models. Although there is clearly scope for refining the model the fit obtained using either of (5.3) or (5.4) is very satisfactory. We hope that our method will allow the possibility of describing accurately the mechanical properties of biological and other polymers with a small number of phenomenological parameters. Moreover, our model can be implemented within standard finite element codes for the analysis of, for example, torsional stresses in DNA loops (e.g. White & Bauer 2004).
We thank Angelo Rosa for providing us with the digitized data used in §6. We are very grateful to Gerhard Holzapfel for his careful reading of the manuscript and for his valuable suggestions. The work of G.S. and I.S. was partially supported by GNFM of Italian INDAM and by PRIN2004 Modelli Matematici per la Dinamica del DNA. RWO is also grateful to GNFM of Italian INDAM for support.