Rationally designed synthetic protein hydrogels with predictable mechanical properties

Designing synthetic protein hydrogels with tailored mechanical properties similar to naturally occurring tissues is an eternal pursuit in tissue engineering and stem cell and cancer research. However, it remains challenging to correlate the mechanical properties of protein hydrogels with the nanomechanics of individual building blocks. Here we use single-molecule force spectroscopy, protein engineering and theoretical modeling to prove that the mechanical properties of protein hydrogels are predictable based on the mechanical hierarchy of the cross-linkers and the load-bearing modules at the molecular level. These findings provide a framework for rationally designing protein hydrogels with independently tunable elasticity, extensibility, toughness and self-healing. Using this principle, we demonstrate the engineering of self-healable muscle-mimicking hydrogels that can significantly dissipate energy through protein unfolding. We expect that this principle can be generalized for the construction of protein hydrogels with customized mechanical properties for biomedical applications.

Then the hydrogels were soaked in PBS, pH 7.4, at the room temperature. After certain time, the hydrogel rings were taken out of PBS buffer, blotted onto tissue paper to remove excess buffer and weighted as W t . The swelling ratio was calculated according to the formula: Swelling ratio (%)＝ (W t -W 0 )/ W 0 ]×100%. Two different samples were measured and the average value was reported. The error bars represent standard deviation.

Supplementary Figure 3 Erosion profiles of 100 mg of Gel 1-4 at room temperature in 100
mM phosphate buffer, pH 7.4. All hydrogels (at 180 mg mL -1 ) were with a surface area of 0.86 cm 2 . The erosion was too slow to be reliably determined using UV spectra if the gel samples were kept still in solution. Therefore, the erosion was measured under mild mechanical shaking.
In the experiment, 100 mg of hydrogel was transferred into a cylindrical glass tube with a flat bottom (1.05cm diameter). The glass tube with the hydrogel was then centrifuged at 1700 g for 10 minutes to completely flat down hydrogel sample to the bottom and smooth the surface of the hydrogel. The hydrogel was allowed to stand overnight. The thin gel film together with the glass tube was then soaked in 5 mL of 100 mM phosphate buffer, pH 7.4, in a scintillation vial. The whole setup was placed on a compact rocker tilting at 50 rpm with amplitude of ±9°, at room temperature. The erosion profiles were determined by measuring the protein absorbance at 280 nm of the supernatant at successive time points using Nano-drop ultraviolet-visible spectrophotometer. Two different samples were measured and the average value was reported. a) Representative gel pictures for Gel-1, Gel-2 and Gel-3 after equilibrium in PBS buffer for 1 day. b) General procedure for the self-healing test. First, the hydrogel ring was completely cut to produce a broken ring. Then, the broken ring was put back to the mold and clamped to move two broken interfaces together for different healing time. After healing, the hydrogel ring became intact and free-standing. The stress-strain curves for the hydrogel rings healed for different time was shown in c, d, and e. The pristine or the self-healed hydrogel was stretched until break at a constant strain rate of 20 mm min -1 . Figure 8 The stretching-relaxation cycles of Gel-2 and Gel-3 after selfhealing for different time periods. All the experiments were at a constant strain rate of ~ 20 mm min -1 . a) The stress-strain curves for Gel-2. The pristine or the self-healed hydrogel was stretched to 140% strain, and then immediately relaxed to 0% strain. b) The stress-strain curves for Gel-3. The pristine or the self-healed hydrogel was stretched to 120% strain, and then immediately relaxed to 0% strain. Because the self-healed Gel-1 can only be stretched to a strain of less than 15%, the stress-strain curves after self-healing were not measured. shape hydrogel upon the same loading as that in experiments with ABAQUS. In the simulation, outer diameter of the ring is 20 mm, the inner diameter of the ring is 16 mm, and the cross-section of the ring is circular. The material is chosen to be neo-hookean with an initial Young's modulus of 130kPa. Due to symmetry, only a quarter of the ring is investigated. The deformed shape and the contour maps of local stresses are provided: (b) Initial shape and deformed shape of a quarter of a soft ring, (c) Contour map of s 1 , (d) Contour map of s 2 , and (e) Contour map of s 3 . The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. Compared with 2 min group: ns, not significant; Compared with 1 d group: NS, not significant. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. Compared with 2 min group: ns, not significant; Compared with 1 d group: NS, not significant. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. ns, not significant. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. ns, not significant, *P < 0.05, **P < 0.01. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. ns, not significant, ***P < 0.001. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. Compared with 2 min group: ns, not significant; Compared with 1 d group: NS, not significant. The hydrogel was stretched until break (failure) at a constant strain rate of 20 mm min -1 . The local slope at a given strain on the stress-strain curve was taken as the Young's modulus at this strain. Student's t-test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant. ns, not significant, *P < 0.05, **P < 0.01.      Table 13 Comparison of some parameters used in the simulation for three different gels, as also listed in Supplementary Tables 10-12 Gel-1 Gel-2 Gel-3

Synthesis of Kir peptide
The Kir peptide with the sequence of CNISYWRESAI was synthesized via the common solid phase peptide synthesis protocol. First, Fmoc-Ile-OH (1 equiv.) and DIPEA (4 equiv.) were dissolved in CH 2 Cl 2 (20 mL g -1 resin). Then chlorotrityl chloride resin (1 equiv.) was added to the solution and the mixture was stirred at room temperature for 1 h. The reaction mixture was filtered and the unreacted resin was capped with 1: 2: mL per gram resin). After the capping procedure, the resin was thoroughly rinsed with CH 2 Cl 2 , DMF, and CH 2 Cl 2 , and then dried in vacuo. The bead-loading was determined to be ~ 0.

Protein expression and purification
The gene encoding polyproteins TIP-1-(GB1) 8  BamHI and BglII. For example, to generate TIP-1-(GB1) 8 -TIP-1, the plasmid for TIP-1 gene was digested with BamHI and KpnI, resulting in an TIP-1 insert with overhanging "sticky ends" whose sequence corresponded to that of the pQE80L vector digested with the same enzymes.
The sticky-ended TIP-1 insert was subsequently ligated into the digested pQE80L vector to form pQE80L-TIP-1. In a similar way, Digestion of (GB1) 8 with BamHI and KpnI resulted in overhanging "sticky ends" whose sequence corresponded to that of the pQE80L-Tip1 vector digested with BglII and KpnI. The sticky-ended (GB1) 8 Figure   11).

AFM tip/Substrate surface modification
The cantilevers for the measurements of the mechanical properties of LBMs were used directly without any modification.
For the measurement of the binding forces for the TIP-1:Kir complex, the cantilever modification procedure was similar as that reported in our previous study [10][11][12][13]  Glass substrates were first immersed into chromic acid for 2 h to remove impurities. After rinsing with deionized water, the glass slides were dried by blowing with nitrogen gas and then placed in a DMSO solution containing silane-PEG-NHS (MW: 5000 Da, 1 mg mL -1 , Nanocs) for 2 h. The glass slides were extensively rinsed with DMSO to remove the unreacted silane-PEG-NHS. Then they were covered by NH 2 -BG solution in DMSO (10 µg mL -1 ) for 2 h so that Snap protein can directly bind to the surface via BG-Snap covalent bond. Finally, Subsequently, the surface was rinsed with deionized water to remove unreacted NH 2 -BG. The glass substrates were used immediately after modification.

Preparation of cys-Xmod-Doc terminated 4-Armed PEG MCL
The protein cys-Xmod-Doc was capped with cysteine at its N-terminus, which can react with

Preparation of Coh-(GB1) 4 -SS-(GB1) 4 -Coh
After the protein Coh-(GB1) 4 -cys has been purified by Co 2+ affinity column, the purified protein was dialyzed against PBS (pH7.4) for 3 d to get Coh-(GB1) 4 -cys-cys-(GB1) 4 -Coh by the formation of the disulfide bond between Coh-(GB1) 4 -cys during dialysis. The dimerization of the protein was very efficient without the need of additional oxidants, as the cys residue after GB1 is highly reactive. Then the resulted product was dialyzed against deionized water (Dialysis membranes MWCO 3, 000; Pierce Co.) for 2 h, repeating for 5 times, to remove salts in the PBS buffer. Finally, the sample was lyophilized to give a white powder of Coh-(GB1) 4 -SS-(GB1) 4 -Coh. The purity was higher than 90% as confirmed by SDS-PAGE.

Hydrogel Preparation
Preparation of the hydrogel was based on the specific protein-peptide interaction between TIP-1 and the Kir peptide or the specific protein-protein interaction between Xmod-Doc and Coh. In figure 3c-e, the hydrogel was stretched until break at a constant strain rate of 20 mm min -1 .
In figure 3f, the hydrogel was stretched to the given length (5% or 10% strain, respectively), and then immediately relaxed to 0% strain. All the experiments were at a constant strain rate of 20 mm min -1 .

Figure 5
In figure 5f, the hydrogel (Gel-4) was stretched until break at a constant strain rate of 20 mm min -1 .
There are totally 5 continuous cycles without any waiting time between each cycle. All the experiments were at a constant strain rate of 20 mm min -1 .
All the experiments were at a constant strain rate of 20 mm min -1 .
In figure 5k, the pristine or the self-healed hydrogel (Gel-4) was stretched until break at a constant strain rate of 20 mm min -1 .
In figure 5l, the pristine or the self-healed hydrogel (Gel-4) was stretched to 160% strain, and then immediately relaxed to 0% strain. All the experiments were at a constant strain rate of 20 mm min -1 .
Supplementary Figure 4 In Supplementary Figure 4, the hydrogel was stretched until break at a constant strain rate of 20 mm min -1 .
There are totally 5 continuous cycles with indicated waiting time between each cycle. All the experiments were at a constant strain rate of 20 mm min -1 .
Immediately after that, the hydrogel was relaxed to the given residual strain following the experimental protocol shown in Supplementary Figure 6a. There are totally 5 continuous cycles without any waiting time between each cycle. All the experiments were at a constant strain rate of 20 mm min -1 .

Supplementary Figure 7
In Supplementary Figure 5c-e, the pristine or the self-healed hydrogel was stretched until break at a constant strain rate of 20 mm min -1 .
Supplementary Figure 8 In Supplementary Figure 6, the pristine or the self-healed hydrogel was stretched to 140% strain, and then immediately relaxed to 0% strain. All the experiments were at a constant strain rate of 20 mm min -1 .

Self-healing Experiment
The broken ring-shaped hydrogel was put back into mold and clamped to bring the two broken interfaces together. PBS buffer was used to cover the ring to avoid dehydration of the ring during the self-healing process. After a preset time, the ring-shaped hydrogel was very carefully take out for further mechanical test. The ring-shaped hydrogel can be partially healed and remain intact when healing time was more than 1 min. The statistics for recovery of the stain and stress of the self-healed hydrogels are summaried in Supplementary Figure 13 and 14, respectively.

Constitutive modeling of hydrogels with folded domains
The synthetic material, schematically shown in Figure 4a, is represented with a volume element of a cube in Fig. 4b. Within this cube, each of 8 protein chains is crosslinked with one arm of linker proteins extending from the cubic center and with that from each corner, essentially similar to the 8-chain model 17 . At the dry state, the cube of RVE is of dimension, l d . Due to solvent absorption or mechanical loading, the represented volume element (RVE) is deformed and can become rectangular within principal axes of stretches. At the current state, the dimensions of RVE become l 1 , l 2 , and l 3 .
To a good approximation, the volume of the RVE at the current state is equal to the sum of the volume of a dry protein network and that of the absorbed water 6 , i.e., where W is the volume per water molecule and M is the number of water molecules.
At the current state, the principal stretches of the cube, denoted as 1 l , 2 l , and 3 l , respectively, are given by where = / . ( is the nominal concentration of water.
When protein chains are subjected to force, F, their originally folded domains can be unfolded.
In single molecule force spectroscopy experiments, the variation of the peak force values is an intrinsic feature because unfolding of folded domains within a polyprotein chain is stochastic.
where α = 1 or Type 1 domain, α = 2 for Type 2 domain, 3* 45 is the unfolding rate at F = 0, ∆ 3 45 is the distance to the transition state on the unfolding pathway, k B is the Boltzmann constant, and T is the absolute temperature.
The unfolded domains can be re-folded with a folding rate, denoted as 3 5 , given by 18,19 0 exp where 3* 5 is the folding rate at F = 0 and ∆ 3 5 is the distance to the transition state on the folding pathway.
The fraction of unfolded domains of each type within a protein chain, L a (a = 1, 2), evolves as The contour length of the protein chain at the current state, denoted as L c , changes due to domains being folded or unfolded and is calculated as where + 5 is the contour length of a chain in the fully folded state, N a (a = 1, 2)is the total number of folded domains of each type within a protein chain, and DL a (a = 1, 2) is the change in contour length when a single folded domain of Type 1 or Type 2 is unfolded.
The force-stretch curve of a protein chain is described by the worm-like chain theory, given by 20 where x is the persistence length of the protein chain and = ( +, − 1) * +, . Two arms of linker proteins in series with a protein chain behaves like a linear elastic spring so that F should also where cl k is the spring constant.
Bonds formed between protein chains and arms of linker proteins may break and re-associate, which has not been considered in the current analysis for simplicity. The failure stress for Gel-1, Gel-2 and Gel-3 is set to be 20 kPa when determing the failure strain in the analysis.
In the theory, the current state of a block is defined by l 1 , l 2 , l 3 , and Λ 3 . It is envisaged that the elastic equilibrium takes place at a much smaller time scale than single stochastic event, including domains being folded or unfolded. This would allow us to decouple the elasticity from the stochastic.
To solve the elastic field, we employ the principle of virtual work. At fixed Λ 3 , we let the RVE at the current state change its dimensions by infinitesimal small amounts, dl 1 , dl 2 , and, dl 3 . The resulting virtual work done by external forces, denoted as P 1 , P 2 , and P 3 respectively, will be Since dl i , i =1, 2, and 3, in Supplementary Eq. 17 are arbitrary and independent variables, we The change in U due to dl i of the RVE is the sum of the change in the elastic energy of protein chains, denoted as U 1 , that in the elastic energy of arms of linker proteins, denoted as U 2 , and that in the energy of mixing water with polymers, denoted as U 3 , i.e., Among them, where N is the number of chains per unit reference volume and +, is the elastic energy of a single protein chain at the current state, given by Since two arms of linker proteins in series behaves as a linear spring, we have ( ) In the theory, we regard that the energy of mixing water with polymers is given by 21,22 where c is a measure of the interaction between polymer and water.
As shown in Figure 3b, circular rings made of gels were loaded by two metallic hooks in experiments. Since gels were very soft and these rings were easily to be bended, very small loads were required to deform rings to a straight configuration in experiments, as shown in the middle of Figure 3b. This straight configuration was then taken as the initial configuration in our experiments. Beyond this point, the hydrogels were regarded to be under uni-axial tension, as seen on the right of Figure 3b. Correspondingly in our theory, we modeled gels with a representative volume element under uniaxial tension along "1" direction with s 2 = s 3 = 0, as illustrated in Supplementary Figure 15a which are set as initial conditions. In Figure 4c-e, the stress is the nominal tensile stress along "1" direction, denoted as s 1 , given by s 1 = s 1 l 2 l 3 /l s2 , and the strain is given by l 2 /l s -1 with the load being applied along "1" direction. For Gel-1 or Gel-3, only one type of domain is considered in the analysis.
The failure of Gel-1, Gel, 2 or Gel-3 in the experiments upon mechanical loadings is due to fracture.
According to the fracture theory 23 , the failure stress, denoted as s f , is given by where E is taken as the initial elastic modulus, G the fracture toughness, and a the crack size.
Considering that the fracture toughness of gels, denoted as G, is due to the rupture of the protein chains lying across the path of the crack 24 , where v is the energy required to rupture a single protein chain and N 2/3 corresponds to area density of protein chains. Note that protein chains are initially subjected to a pre-stretch at the free swelling state of gels. The final rupture of a single protein chain is due to bond breaking of one arm of linker proteins that is in series with the protein chain instead. Due to the nature of multiple bonds existing within gels, the apparent bond breaking force of a crosslinker should also depend on its rebinding kinetics. In calculating v, a protein chain is subjected to a constant stretching velocity, 10 nm s -1 and the apparent bond breaking force is set to be 20 pN, based on the loading rate dependent experiments for the rupture forces shown in Extended Data Figure 1. v is found to be 92.3 zJ for Gel-1, 464 zJ for Gel-2, and 977 zJ for Gel-3, respectively. With Supplementary Eq. 26-27, the ratio among failure stresses for Gel-1, Gel-2, and Gel-3 is calculated to be 1.4:1:2, which is close to the experimental value, 1:1:2, shown in Figs. 3c-d. The crack size in gels is also estimated, which is on the order of 0.5 µm.
The RVE used in our theory was essentially similar to the 8-chain model 17 , which was a classical network model to describe mechanical properties of polymeric materials. The 8-chain model 17 was mathematically simple but also able to capture various deformation modes, including uniaxial tension, pure shear, and biaxial tension. By now, the 8-chain model 17 has also been incorporated into various theories to describe gel deformation 25,26 . Other polymeric models, such as 3-chain model 27 , 4-chain model 28,29 , or full network model 30 may also model our system well and the similar trends can be expected. The 8-chain model we used can capture shear deformation well, as documented in literature 31 .
Among parameters adopted in our simulation as listed in Supplementary Tables 10-12, x, W, * +H , k cl , g, c, µ were the same for Gel-1, Gel-2, and Gel-3. The unfolding rate, the unfolding distance, the refolding rate, or the refolding distance for different domains within LBMs for Gel-1, Gel-2, Tables (10)(11)(12), are also given in Supplementary Table 13 for better comparison. It is expected that all these parameters would affect simulated rigidity, toughness and extensibility in our theory.

The comparison between theory and experiments
The comparison of Young's modulus, failure strain, and fracture toughness between theory and experiment was provided in Supplementary Tables 14-16, respectively. Note that the fracture toughness from experiments was obtained through Eq. (S26) with measured E and s f from experiments and also an assumption of a crack size of 0.5 µm. As seen from Supplementary Tables 14-16, there exist mismatches between experiment and theory.
In our theory, it was regarded that proteins gels form a crosslinked three-dimensional network with its mechanical properties described by RVE shown in Fig. 3b. However, in reality, the crosslinked protein network within fabricated gels might not be so perfect as that illustrated in Figure 1a. For example, some arms of CLs might not be able to form crosslinks with LBMs and some LBMs might have just hanged to the network from one of its end or both ends might be separated from the network. Due to these defects, protein length in the main network may be inhomogeneous and unfolding or re-folding dynamics of folded domains with LBMs would be affected, which are expected to affect mechanical properties of fabricated gels in turn. In addition, the loading condition was assumed to be uniaxial tension in the simulation, which is slightly different from that in our experiments. These might be part of the reasons why the shapes of hysteresis obtained from experiments look different from our theoretical prediction, as comparing Figure 3 with Figure   4.
It should be pointed out that the crack size, denoted as a in Eq. (S26), was not determined. We therefore suggest that the discrepancy of facture toughness between theory and experiment might be due to the different crack sizes existing within three different gels. Nevertheless, defects in the crosslinked protein network, voids, and other damages, as well as their evolution under loads might also have contributed this discrepancy.
To improve our theory further, inhomogeneity of protein network of fabricated gels may be captured by varying the length, number of domains, or unfolding or re-folding dynamics of domains, etc. of 8 chains within a RVE. To better predict fracture toughness, it is desirable to find out initial defects within gels in experiments. These defects should then be incorporated into a fracture theory, which also considers how these defects evolve. In addition, the breaking and reformation of bonds between LBMs and CLs have not considered in our theory yet, which should be included in predicting the healing process of gels in our future work. It would also be great to simulate the load-displacement curve of hoop samples, as used in the experiments, by incorporating the constitutive relation developed in our work into a program for finite element method, such as ABAQUS.

Statistical analysis
Igor Pro and Microsoft Excel were used to plot and analyze all the graphs. Student's t test was used for statistical analyses. When P value is below 0.05, it is considered being statistically significant.