Transient viscoelasticity study of tobacco mosaic virus/Ba2+ superlattice

Recently, we reported a new method to synthesize the rod-like tobacco mosaic virus (TMV) superlattice. To explore its potentials in nanolattice templating and tissue scaffolding, this work focused the viscoelasticity of the superlattice with a novel transient method via atomic force microscopy (AFM). For measuring viscoelasticity, in contrast to previous methods that assessed the oscillating response, the method proposed in this work enabled us to determine the transient response (creep or relaxation) of micro/nanobiomaterials. The mathematical model and numerical process were elaborated to extract the viscoelastic properties from the indentation data. The adhesion between the AFM tip and the sample was included in the indentation model. Through the functional equation method, the elastic solution for the indentation model was extended to the viscoelastic solution so that the time dependent force vs. displacement relation could be attained. To simplify the solving of the differential equation, a standard solid model was modified to obtain the elastic and viscoelastic components of the sample. The viscoelastic responses with different mechanical stimuli and the dynamic properties were also investigated.


Background
The recognition of tobacco mosaic virus (TMV) since the end of nineteenth century [1] has sparked innumerable research towards its potential applications in biomedicine [2,3] and biotemplates for novel nanomaterial syntheses [4,5]. A TMV is composed of a single-strand RNA that is coated with 2,130 protein molecules, forming a special tubular structure with a length of 300 nm, an inner diameter of 4 nm, and an outer diameter of 18 nm [6]. The TMVs observed under a microscope can reach several tens of microns in length due to its unique feature of head-to-tail self-assembly [7]. Practically useful properties of the TMVs include the ease of culture and broad range of thermal stability [8]. Biochemical studies have shown that the TMV mutant can function as extracellular matrix proteins, which guide the cell adhesion and spreading [8]. It has also been confirmed that stem cell differentiation can be enhanced by both native and chemically modified TMV through regulating the gene's expression [9][10][11]. Moreover, TMV can be electrospun with polyvinyl alcohol (PVA) into continuous TMV/PVA composite nanofiber to form a biodegradable nonwoven fibrous mat as an extracellular matrix mimetic [12].
Very recently, we have reported that the newly synthesized hexagonally packed TMV/Ba 2+ superlattice material can be formed in aqueous solution [13,14]. Figure 1 shows the schematic of the superlattice formation by hexagonal packing of TMVs, triggered by Ba ions, and the images observed from field emission scanning electron microscopy (FESEM) and atomic force microscopy (AFM). The sample we used for this experiment was tens of microns in length, 2~3 microns in width (from FESEM), and several hundred nanometers in height (from AFM height image). It is known that the superlattice exhibits physical and mechanical properties that differ significantly from its constituent materials [15][16][17][18][19][20]. The study on the viscoelastic properties of the TMVderived nanostructured materials is still lacking despite the availability of the elastic property of the TMV and TMV-based nanotube composites [7]. The viscoelasticity of micro/nanobioarchitecture significantly affects the tissue regeneration [21] and repair [22], cell growth and aging [23], and human stem cell differentiation [24] as well as the appropriate biological functions of the membranes within a specific nanoenvironment [25]; in particular, the viscoelasticity of some viruses plays key roles in the capsid expansion for releasing nucleic acid and modifying protein cages for vaccine delivery purposes [26]. Specifically, for TMV superlattice, its nanotube structure makes it a perfect biotemplate for synthesizing nanolattices that have been confirmed to possess extraordinary mechanical features with ultralow density [27,28]. Considering the biochemical functions of the TMV, its superlattice is an excellent candidate for bone scaffolding where the time-dependent mechanical properties become determinant [29], and research on scaffolding materials remains a hotspot [30]. Apart from contributing to the application of TMV superlattice, this work also pioneered in the viscoelasticity study of virus and virus-based materials. By far, most literature on viral viscoelasticity has been focused on the dynamic properties of virus suspensions or solutions [31][32][33][34]. One of the rare viscoelasticity studies on individual virus particle is the qualitative characterization of the viscoelasticity of the cowpea chlorotic mottle virus [26] using quartz crystal microbalance with dissipation technique, which presents only the relative rigidity between two samples. To date, little literature is available on the quantitative study of the viscoelasticity of individual virus/virus-based particles. Considering the potential uses of TMV/Ba 2+ superlattice, its viscoelastic properties and responses under different mechanical stimuli need to be investigated.
A number of techniques for measuring the viscoelasticity of macro-scale materials have been used. A comprehensive review of those methods can be found in the literature [35] that addresses the principles of viscoelasticity and experimental setup for time-and frequency-domain measurements. When the sample under investigation is in micro or even nanometer scale, however, the viscoelastic measurements become much more complicated. In dynamic methods, shear modulation spectroscopy [36] and magnetic bead manipulation [37] are two common methodologies to obtain the micro/nanoviscoelastic properties. To improve the measurement accuracy, efforts have been made to assess the viscoelasticity of micro/nanomaterials using contact-resonance AFM [38][39][40][41]. The adhesion between the AFM probe tip and sample, however, is usually neglected. Furthermore, in order for the dynamic method to obtain a sinusoidal stress response, the applied strain amplitude must be kept reasonably small to avoid chaotic stress response and transient changes in material properties [42]. In addition, the dynamic properties are frequency dependent, which is time consuming to map the viscoelasticity over a wide range of frequencies. An alternative way to measure the viscoelastic response of a material is the transient method. Transient indentation with an indenter was developed based on the functional equation methods [43], where the loading or traveling histories of the indenter need to be precisely programmed.
In this study, the viscoelastic properties of the TMV/ Ba 2+ superlattice were investigated using AFM-based nanoindentation. AFM has the precision in both force sensing and displacement sensing, although it lacks the programing capability in continuous control of force and displacement. To realize the transient indentation in AFM, we introduced a novel experimental method. Viscoelastic nanoindentation theories were then developed based on the functional equation method [44]. The adhesion between the AFM tip and the sample, which significantly affected the determination of the viscoelastic properties [45], was included in the indentation model [20]. The viscoelastic responses of the sample with respect to different mechanical stimuli, including stress relaxation and strain creep, were further studied. The transition from transient properties to dynamic properties was also addressed.

Methods
The TMV/Ba 2+ superlattice solution was obtained from the mixture of the TMV and BaCl 2 solution (molar ratio of Ba 2+ /TMV = 9.2 × 10 4 :1) as stated in the reference [13]. It was further diluted with deionized water (volume ratio 1:1). A 10-μL drop of the diluted solution on a silicon wafer was spun at 800 rpm for 10 s to form a mono-layer dispersion of the sample. The sample was dried for 30 min under ambient conditions (40% R.H., 21°C) for AFM (Dimension 3100, Bruker, Santa Barbara, CA, USA) observation and subsequent indentation tests.
The sample was observed with FESEM and AFM. The indentation was performed using the AFM nanoindentation mode (AFM probe type: Tap150-G, NanoAndMore USA, Lady's Island, SC, USA). The geometry of the cantilever was precisely measured using FESEM (S-4700, Hitachi, Troy, MI, USA), with a length of 125 μm, width of 25 μm, and thickness of 2.1 μm. To accurately measure the tip radius, the tip was scanned on the standard AFM tip characterizer (SOCS/W2, Bruker) and the scanned data was curve fitted using PSI-Plot (Poly Software International, Orangetown, NY, USA). The tip radius calculated to be 12 nm. For a typical indentation test, the tip was pressed onto the top surface of the sample until a predefined force of~100 nN. The cantilever end remained unchanged in position during the controlled delay time. A series of indentations of the same predefined indentation force and different delay times were performed to track the viscoelastic responses. A 10-min time interval of the two consecutive indentations was set for the sample to fully recover prior to the next indentation. The sample drift was minimized by turning off the light bulb in the AFM controller during scanning to keep the AFM chamber temperature constant and by shrinking the scan area gradually down to 1 nm × 1 nm on the top surface of the sample to rid the scanner piezo of the hysteresis effect.

Mathematical formulation
Derived from the functional equation method and the standard solid model (shown in the ' Appendix'), the differential equation governing the contact behavior of viscoelastic bodies can be obtained as where F(t) is the contact force history, δ(t) is the indentation depth history, R is the nominal radius of the two contact spheres, w is the adhesive energy density, A i and B i (i = 0, 1, 2) are the parameters determined by the mechanical properties of two contact bodies, and the calculation of all these parameters can be found in the ' Appendix.' The elastic moduli E 1 and E 2 and viscosity η in Figure 2 are implicitly included in the above differential equation.
To determine E 1 , E 2 , and η, besides experimental data for t and F, the function of the force history F(t) is also required. The experimental data of t and F can be obtained as indicated in Figure 3. The force relaxation can be found in Figure 3a where the force decrease between the right ends of extension and retraction curves. By mapping the force decrease at different delay times as shown using the red asterisks in Figure 3b, the force relaxation curve can be obtained, which decreases from 104 to 40 nN. The function of F(t) can be obtained from Equation (1). Not only is Equation (1) applicable for the standard solid model in Figure 2(a) where it is derived from, but also it can be used for the modified standard solid model in Figure 2(b) where the elastic component of E 1 is replaced by two elastic components in series. With this modification, the deflection of the cantilever can be incorporated into the deformation of the imaginary sample which is represented by the modified standard solid model where the elastic component of E 1c in Figure 2(b) denotes the cantilever and the rest components denote the TMV/Ba 2+ superlattice.
During each indentation, the vertical distance between the substrate and the end of the cantilever remains constant. Therefore, as the sample deformation or the indentation depth increases, the corresponding cantilever deflection Δd or the normal indentation force decreases. During this process, the force on the system decreases while the sample deformation δ increases to compensate the decreased cantilever deflection. Therefore, the change of the cantilever deflection is equal to change of the sample deformation during indentation, as is shown in Figure 4. As such, δ in Equation (1) represents the relative approach between the cantilever end and the substrate, which incorporates the deformation of both the sample and the cantilever.
To be clearer, δ is substituted by D which represents the combined deformation. The relative approach, D, can be written as where H(t) is the Heaviside unit step function and D 0 is the relative approach between the substrate and the end of the cantilever. Thus, Equation (1) can be rewritten as Applying Laplace transform, it yields where a function with ' ∧ ' denotes Laplace-transformed function in s domain. Performing inverse Laplace transform, the viscoelastic equation of AFM-based indentation becomes  where

Solution to AFM-based indentation equation
It is observed from Figure 3 that the initial indentation force at t = 0 was measured to be 104.21 nN, then the force started to decrease and then remained constant at 38 nN after~5,000 ms. The force decrease shown as red asterisks in Figure 3b fits qualatitatively well with the exponential function of Equation (5). E 1 , E 2 , and η, corresponding to the mechanical property parameters in Figure 2(a), can then be determined by fitting Equation (5) with the experimental data. From the indentation data, D 0 is obtained to be 78.457 nm. The pull-off force, 2πwR, calculated by averaging the pull-off forces of multiple indentations on the sample, is 16 nN. In comparison with the radius of the AFM tip, the surface of the sample can be treated as a flat plane. Hence, the nominal radius R = R tip = 12 nm.
By invoking the force values at t = 0, t = ∞, and any intermediate point into Equation (5), the elasticity and viscosity components can be determined to be E 1 = 32.0 MPa, E 2 = 21.3 MPa, and η = 12.4 GPa ms. The coefficient of determination R 2 of the viscoelastic equation and the experimental data is~0.9639.
Since the stress relaxation process is achieved by modeling a combination of the cantilever and the sample, the viscoelasticity of the sample can be obtained by subtracting the component of the cantilever from the results. The cantilever, acting as a spring, is in series with the sample, represented by a standard solid model. The schematic of the series organization is shown in Figure 2 (b). Thus the component of E 1 comprises of E 1s representing the elastic part from the sample and E 1c representing the elastic part from the cantilever. To clarify the sources of the components in the modified standard solid model, E 2 , v 2 , and η in Figure 2(a) are now respectively denoted by E 2s , v 2s , and η s in Figure 2(b), where the subscript 's' denotes the sample.
At the onset of indentation, only the spring with elastic modulus of E 1 takes the instantaneous step load; therefore, the elastic modulus of E 1s can be determined from the experimental data of zero-duration indentation. Applying the DMT model [46] with the force-displacement relationship of the cantilever, we can obtain the elastic equation of AFM-based indentation where k is the spring constant of the cantilever, which is 5 nN/nm based on Sader's method [47] to calibrate k, δ cantilever is the cantilever deflection, and δ is recorded directly as the Z-piezo displacement by AFM. The elastic modulus of E 1s can be calculated by fitting the DMT-model-based indentation equation with experimental data as shown in Figure 5. For simplicity, modification was done to the indentation equation and the experimental data, whose details can be found in reference [20]. The fitted elastic modulus of E 1s is~2.14 GPa with a coefficient of determination of 0.9948.

Results and discussion
Based on the solution obtained, the viscoelastic equation of AFM-based indentation for TMV/Ba 2+ superlattice is written as The force decrease curve is shown in Figure 3b with the experimental data.
Specifically, for the TMV/Ba 2+ superlattice whose viscoelastic behavior is simulated by a standard solid model, the differential equation governs its stress-strain behavior and becomes where E 1s = 3 GPa, E 2s = 21.3 MPa, and η s = 12.4GPa ms.
In the standard solid model, the initial experimental data point is determined by the instantaneous elastic modulus E 1s . For the indentation that is held for over 5,000 ms, the indentation force becomes steady at~38 nN, when the force exerts on the two springs in series. In contrast to E 1s , E 2s is much smaller, as can be seen from the significant force decrease of from~104 to~38 nN. The tip traveled down 13.2 nm from the beginning of indentation. It is noted that for our indentation test, the ratio of the maximum indentation depth to the sample diameter is less than 10% [48,49]; the substrate effect to the elastic modulus calculation is neglected.
From the determined viscoelastic model, the mechanical response of the superlattice under a variety of mechanical loads can be predicted. Several simulation results were included as follows.
When the TMV/Ba 2+ superlattice sample undergoes a uniformly constant tensile/compressive strain, the stress relaxation can be obtained from the standard solid model as below where ε 0 is the constantly applied strain. When the sample undergoes a uniformly constant tensile/compressive stress, the strain creep can then be obtained as where σ 0 is the constantly applied stress. The stress relaxation vs. applied strains and the strain creep vs. applied stresses are shown in Figure 6a,b, respectively. In Figure 6a, the stress reduces to a steady state after~2 s when the applied strain is~10%. In Figure 7b, strain increases to a steady value after~5 s when the applied stress is~1 GPa.
When the sample is indented with a spherical indenter, the indentation depth history can be analytically obtained when a step force is applied. Similar to the procedures above where the force history of Equation (5) is obtained, a step force function is used as input, and the creep indentation depth history function can be derived as ð12Þ where F 0 is the step force, A 0 ¼ 4G 1s G 1s þG 2s þ 3K 1s G 2s Figure 5 Indentation force data as a function of Z-piezo displacement, a comparison of experimental measurement and fitted results.
The indentation force history has been obtained in Equation (5), where the elastic shear modulus G 1 as a combined elastic response of two springs shown in Figure 2(b) should be replaced by G 1s of one spring only. Then, the simulated curves for the two situations can be found in Figures 6c,d. It is concluded that the creep depth variation under different forces gets larger through creep while the indentation force variation under different depths gets smaller through relaxation. Particularly, in Figure 6d, the force finally decreases to negative values, which represent attractive forces. The attraction cannot be found when G 1s and G 2s are very small. This phenomenon can be interpreted by the conformability of materials determined by the elastic modulus. When G 1s and G 2s get smaller, the materials are more conformable. Accordingly, in the final equilibrium state, the materials around the indenter tend to be more deformable to enclose the spherical indenter. This will result in a smaller attraction.
In addition, the example of shear dynamic experiment is simulated to obtain the storage and loss moduli of  TMV/Ba 2+ superlattice. The storage and loss shear moduli are calculated by [42] where G′ and G″ are storage and loss moduli, respectively, ω is the angular velocity which is related to the frequency of the dynamic system, and G s t ð Þ ¼ G 1s þ G 2s e −G 2s t=η is the shear stress relaxation modulus, determined by the ratio of shear stress and constant shear strain.
Based on the relation between the transient and dynamic viscoelastic parameters in Equations (13) and (14), the storage and loss shear moduli are finally determined to be where G 2s = E 2s / 2(1 + v 2s ). Figure 7 shows the curves of storage and loss shear moduli vs. the angular velocity. The storage shear modulus, G′, increases with the increase of angular velocity, while the increasing rate of G′ decreases and the angular velocity of~2 rad/s is where the increasing rate changes most drastically. However, the loss shear modulus, G″, first increases and then decreases reaching the maximum value,~3.9 MPa, at the angular velocity of~0.7 rad/s. The storage and loss moduli in other cases as uniform tensile, compressive, and indentation experiments can also be obtained.

Conclusions
This paper presented a novel method to characterize the viscoelasticity of TMV/Ba 2+ superlattice with the AFMbased transient indentation. In comparison with previous AFM-based dynamic methods for viscoelasticity measurement, the proposed experimental protocol is able to extract the viscosity and elasticity of the sample. Furthermore, the adhesion effect between the AFM tip and the sample was included in the indentation model. The elastic moduli and viscosity of TMV superlattice were determined to be E 1s = 2.14 GPa, E 2s = 21.3 MPa, and η s = 12.4 GPa•ms. From the characterized viscoelastic parameters, it can be concluded that the TMV/Ba 2+ superlattice was quite rigid at the initial contact and then experienced a large deformation under a constant pressure. Finally, the simulation of the mechanical behavior of TMV/Ba 2+ superlattice under various loading cases, including uniform tension/compression and nanoindentation, were conducted to predict the mechanical response of sample under different loadings. The storage and loss shear moduli were also demonstrated to extend the applicability of the proposed method. With the characterized viscoelastic properties of TMV superlattice, we are now able to predict the process of tissue regeneration around the superlattice where the time-dependent mechanical properties of scaffold interact with the growth of tissue.

Modeling of adhesive contact of viscoelastic bodies
The functional equation method was employed to develop a contact mechanics model for indenting a viscoelastic material with adhesion. A modified standard solid model was used to extract the viscous and elastic parameters of the sample.
Several adhesive contact models are available, such as Johnson-Kendall-Roberts (JKR) model [50], Derjaguin-Muller-Toporov (DMT) model [46], etc. [51][52][53]. Detailed comparisons can be found in reference [54]. As the DMT model results in a simpler differential equation, it was used in this study for the simulation to solve the indentation on an elastic body with adhesion.
For the DMT model [46], the relation between the indentation force F and relative approach δ, shown in Figure 8, can be expressed as where R is the nominal radius of the two contact spheres of R 1 and R 2 , given by R = R 1 R 2 /(R 1 + R 2 ); the adhesive energy density w is obtained from the pull-off force F c , where F c = 3πwR/2; and the reduced elastic modulus E * is obtained from the elastic modulus E s and Poisson's ratio ν s of the sample by with the assumption that the elastic modulus of the tip is much larger than that of the sample.
In Equation (A.1), E * , which governs the contact deformation behavior, is decided by the sample's mechanical properties. In the functional equation method [43], E * needs to be replaced by its equivalence in the viscoelastic system, so that the contact deformation behavior can be governed by the viscoelastic properties. To achieve it, the elastic/viscoelastic constitutive equations are needed.
As a premise of the functional equation method, quasi-static condition is assumed so that the inertial forces of deformation can be neglected [43,44]. The general constitutive equations for a linear viscoelastic/elastic system in Cartesian coordinate configuration can be written as ðA:2Þ where s ij , e ij , σ kk , and ε kk are the deviatoric stress, strain, mean stress, and strain, respectively. The linear operators P d , Q d , P m , and Q m can be expressed in the form of where i (i = 0, 1, 2,…) is determined by the viscoelastic model to be selected, t is time, and p d i , q d i , p m i , and q m i are the components related to the materials property constants, such as elastic modulus and Poisson's ratio etc.
For a pure elastic system, the four linear operators are reduced to which, according to the elastic stress-strain relations, are correlated as where G and K are the shear modulus and bulk modulus, respectively.
Combining Equation (A.6) with the reduced elastic modulus can be expressed by the elastic linear operators as To evolve the elastic solution into a viscoelastic solution, the linear operators in the viscoelastic system need to be determined. To this end, the standard solid model, shown in Figure 2(a), was used to simulate the viscoelastic behavior of the sample, since both the instantaneous and retarded elastic responses can be reflected in this model, which well describes the mechanical response of most viscoelastic bodies.
It is customary to assume that the volumetric response under the hydrostatic stress is elastic deformation; thus, it is uniquely determined by the spring in series [55]. Hence, the four linear operators for the standard solid model can be expressed as Þ , E 1 , E 2 , v 1 , and v 2 are the elastic modulus and Poisson's ratio of the two elastic components, respectively, shown in Figure 2.