The DSB coarse-grained model

The Dynamic Structure-Based model has been introduced in ref. [9] and here we just give a brief summary. We label the interacting pseudoatoms by indices i and j. Their positions are $$ {\overrightarrow{r}}_i$$ and $$ {\overrightarrow{r}}_j$$, and $$ r_{i,j}=|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j|$$ is the distance between them. The mass of each residue is set to the average amino acid mass m. The positions evolve according to the Langevin equations of motion with the time unit, τ, of ≈ 1 ns, the damping coefficient γ = 2m/τ (corresponding to overdamped dynamics [24]), and thermal white noise $$ \overrightarrow{{\Gamma }_i}$$ with variance σ² = 2γk_BT:

We solve the equations of motion by using the 5th order predictor-corrector algorithm [25]. Energy unit ϵ is of order 1.5 kcal/mol [26, 27]. The room temperature is then about 0.3–0.35 ϵ/k_B, where k_B is the Boltzmann constant. Here, the simulations are run at 0.3 ϵ/k_B. The results are presented, approximately, in the metric units.

The excluded volume of the residues (whether involved in a contact or not) is ensured by the L-J potential that is cut at r_o = 0.5 nm (so that V_r(r_i,j ≥ r_o) = 0):

The distinction between the bb, bs and ss contacts in our one-bead-per-residue model is accomplished in the following way. When two beads come sufficiently close to each other, we compute two vectors based on positions of three consecutive beads (i-1,i,i+1): $$ -{\overrightarrow{n}}_i$$ is the negative normal vector, approximately pointing in the sidechain direction, and $$ {\overrightarrow{h}}_i$$ is the binormal vector that points in the direction of a possible backbone hydrogen bond. The relative directions of these vectors are checked, and if, for instance, $$ {\overrightarrow{h}}_i$$ and $$ {\overrightarrow{h}}_j$$ point towards each other, a bb contact is made. Each type of amino acid has a specific number of contacts it can make with its virtual “backbone” and “sidechain”. Each ss contact making a pair of amino acids corresponds to different distances at the minimum of the L-J potential. On the other hand, that minimum is 0.5 nm for the bb contacts and 0.68 nm for the bs contacts. Independent of the source of the contact, the interaction is modeled by the potential of the same form:

Simulation protocol

We performed simulations of 5 different systems specified in Table 1. The exact protein composition of all simulated systems is presented and described in detail in Section 1 of the S1 Text.

While cereal storage proteins are still in the grain, they are usually located in the protein bodies or protein storage vacuoles [38]. After grain maturation, these deposits can fuse and form a continous network [1]. This is reflected in our simulation protocol: in the beginning, the proteins are in a low-density box and have time to attain their shapes separately [39]. Then, as the box dimensions get smaller, different chains connect together, forming one big cluster. Unfortunately we cannot reproduce drying and hydratation in an implicit-solvent model, but the final density of the system reflects the density of real gluten samples [40]. Mixing and kneading of the flour and the resulting dough is represented by the periodic deformations of the box, that can also be compared to the experimental deformations of (much larger) gluten samples in a rheometer [41].

The first step of the simulations is generating the chains as self-avoiding-random-walks positioned randomly in a box with density ρ = 0.1 res/nm³ (residues per cubic nanometer). If an initial random walk configuration happens to cross the boundary of the box, the box is enlarged to keep all chains inside, resulting in a slighly lower initial density. The density of 0.1 res/nm³ is more than 30 times lower than the density of gluten, which should ensure that proteins are in monomeric or dimeric form during the initial equilibration (see inset 1 in Fig 1). Boundary conditions along the X and Y axes are periodic, whereas the walls along the Z direction are repulsive. They are separated by the distance of s. The residues are repelled from the wall by the repulsive part of the L-J potential with the depth 4ϵ.

Fig 1

Time dependence of the distance s separating the two opposite box walls in the Z direction.

The blue curve shows the force exerted on the system by those two walls (summed over all the residues and averaged over 100 ns time interval). The data is for gluten (defined in Section 1 of the S1 Text), with the density of 3.5 nm⁻³. The period of oscillations is 40 ms. The six insets correspond to the snapshots of the system taken at different stages of the simulation (each bead represents one amino acid, the protein chains are shown in different colors). The same coordinate system (top left) is used for all 6 insets.

After equilibrating the system for 125 μs the box dimensions are reduced from the initial size with the speed of 2 mm/s, after the density ρ₀ (usually corresponding to the real density of gluten [40] ρ₀ ≈ 3.5 res/nm³) is reached (see Section 3 of the S1 Text for more information about the density). At this stage, the distance between the repulsive walls has a value that is denoted by s₀. The system is then equilibrated for another 150 μs.

After the second equilibration stage (depicted in inset 2 in Fig 1) the walls along the Z direction start to attract the proteins. When the distance between any residue and the wall gets closer than d_min = 0.5 nm, the attractive part of the L-J potential with depth of 4ϵ is turned on adiabatically (in the same manner as for disulfide bonds). Then the interaction between the wall and a residue attached to it is described by V_wall = 4V^LJ(r_i,w), where r_i,w is the distance between the ith pseudoatom and an ith interaction center on the wall. The X and Y coordinates of the ith interaction center are fixed at the values of the first bead attachment (when it came to the wall closer than d_min = 0.5 nm, which is the distance corresponding to the minimum of the V_wall potential). The Z coordinate of an interaction center is 0 on one rigid wall and s on the opposite wall. If a wall moves so that s ≠ s₀, the interaction centers move together with the wall. If the bead departs away from the interaction center by more than 2 nm, it is considered detached and the interaction center disappears. The bead may reattach at another place later. For d < d_min the repulsive part of the L-J potential is always turned on to ensure that no residue can pass through the wall. The force acting on the walls is a sum of forces that individual residues exert on both walls. It is averaged over 100 ns interval to filter the thermal noise out. An alternative way to introduce interactions with the walls (involving the fcc lattice of psudoatoms) is discussed in Section 4 of the S1 Text.

After the wall attraction is enabled, we give another 150 μs for the equilibration of the system (inset 3 in Fig 1). After this third equilibration stage, when the attraction is already turned on, the walls in the Z direction are moved to the distance of s = s₀ − A, where A is the amplitude of the oscillations in the normal direction (we set A = 1 nm). After being moved to the minimum (shown in inset 4 in Fig 1), the wall position changes periodically: s(t) = s₀ − A cos(ωt). The wall at Z = 0 stays fixed. In the experiments, ω ∼ 1 Hz [2], which is not possible to achieve in our simulations. Our default period is 40 μs, corresponding to f ≈ 25 kHz. The maximum displacement is shown in inset 5 in Fig 1. For shearing simulations, the walls are also displaced with the same amplitude, A, but the oscillations take place in the X direction and are described by s′(t) = −A cos(ωt) (where s′ = 0 corresponds to one rigid wall being directly above the other). See S1 Fig for an illustration.

After oscillations in the normal direction (or after additional 100 μs of equilibration if no oscillations were made) the system is equilibrated again for 75 μs and then the walls in the Z direction extend gradually to 2s₀ so that the total work required for such an extension can be measured, as well as the distance and force required to break the system apart (recreating conditions for measuring uniaxial extension [41]). Inset 6 in Fig 1 depicts this elongation. The speed of elongation is 0.5 mm/s, which is a typical pulling speed employed in AFM stretching experiments in coarse-grained simulations [26, 42]. The maximum speed during oscillations (v_max = Aω) is even lower.

Properties of gluten, maize and rice

Based on refs. [3, 8, 11, 43–49] we chose proteins that were representative for each of the five systems we study: gliadins, glutenins, gluten and storage proteins from maize and rice. We then used the UniProt database [35] to select the corresponding sequences. The resulting compositions are described in Section 1 of the S1 Text. Table 1 summarizes the key properties of each of the systems studied. A large number of cysteines per chain should increase the propensity to form inter-chain disulfide bonds. On the other hand, longer chains should make entanglements easier. Finally, the “stickiness” of the system is reflected in the coordination number z, which is defined as the mean number of contacts made by each residue (both “dynamic” and “static” combined).

We performed three types of simulations: a) without any oscillations; b) with oscillations with deformation in the normal (Z) direction; c) with oscillations with deformation in the shear (X) direction. All of the properties are calculated in the stage after oscillations (in the absence of oscillations—after the equivalent time is passed). The parameter z does not seem to be affected by the oscillations.

However, the values of 5 other quantities depend on the prior history of the system. Those quantities are: the number of entanglements l_k, the maximum force F_max, the mechanical work W_max required to stretch the system, the number of the interchain contacts n_inter and the RMSF (root mean square fluctuation) as averaged over all residues and determined as a deviation from a mean location of the residue.

It should be noted that each of the systems studied is simulated in a somewhat different box (to achieve the same density for different numbers of residues) and thus comparison of the absolute values of these quantities is not very meaningful. Instead, we analyze the ratios of the five quantities to their average values in the absence of oscillations (denoted by the tildas). Fig 2 shows these ratios for the case of the shearing oscillations whereas S3 Fig for the uniaxial oscillations.

Fig 2

The ratio of 5 different properties calculated after 5 shearing oscillations (in the nominator), and without any oscillations (in the denominator, marked by ~).

The properties are the average number of entanglements l_k, the maximum force F_max during the elongation in the last simulation stage and the maximum work W_max required to elongate the system in that stage, the number of inter-chain contacts n_inter and RMSF (root mean square fluctuation) averaged over all residues. The ratio of 1 is marked by a broken line. Each color corresponds to a different system, as displayed at the top, and defined in Section 1 of the S1 Text. The system density is 3.5 nm⁻³.

The first of the quantities studied is the total number of entanglements. We have determined l_k by using the Z1 algorithm [17]. An entanglement occurs if a path connecting the first and last residue in a chain cannot be shortened to a straight line without crossing a path from another chain [18, 19]. After minimizing the lengths of all such paths with the constraint that they cannot cross one another, each kink, where two paths intertwine, counts as an entanglement [20]. An example of an entanglement is shown in inset 6 in Fig 1, where the end of the grey chain (a low molecular weight glutenin) crosses a loop made by a blue chain (a high molecular weight glutenin). Due to the low number of the chains, each snapshot of the simulation usually has only a few entanglements. We averaged l_k over all snapshots made after oscillations.

We use a deterministic random number generator, which produces the same output when provided with same initial random number, called “random seed” [50]). The simulations for each system were repeated at least three times with a different random seed. The same set of random seeds was used for the shearing, pulling and equilibrium simulations. The result shown in Fig 2 uses the mean values from all simulations, but the error of the mean is quite high. Nevertheless, a clear increase is observed for gluten: the oscillating deformations cause more entanglements. We did not include self-entanglements (knots) in Fig 2, but they are also present in gluten: some long proteins can form very deep knots (see S5 Fig), that can last for over 100 μs.

The differences between the studied systems are most pronounced in the case of F_max, the maximum force required to stretch the system in the last stage of the simulation (illustrated in inset 6 in Fig 1). The blue curve in Fig 1 shows how the force rises and then drops near the end due to the rapture. However, different random seeds may lead to different force curves for the same system. Fig 3 shows three stretching force curves for three different simulations of gluten. Each curve has a number of force peaks. In order to determine the characteristic largest force, F_max, we take 5 largest force points and calculate their mean value. These mean F_max are also indicated in Fig 3 as dotted lines surrounded by strips that represent their standard deviation. The F_max used in Fig 2 is a weighted average of those values (from simulations with different random seeds), where the weights are calculated from the standard deviation of the mean. Large glutenins experience the largest, 5-fold increase in F_max after oscillations, while for gliadins there is even a very slight decrease. The rise in F_max for gluten is substantially bigger than for corn and rice proteins.

Fig 3

The force exerted on gluten with the density of 3.5 nm⁻³ (top panel) during elongation in the last stage of the simulation and the resulting work (bottom panel) as a function of the distance between the rigid walls.

The three trajectories are shown in different colors. The dotted lines correspond to F_max and W_max for a given trajectory. For F_max the strips around the dotted lines represent the uncertainty resulting from averaging F_max over 5 points with the biggest force. The oscillation period is 40 μs. For F, the thin lines represent raw simulation data and the thick lines are smoothed.

The last stage of the simulations includes moving the Z walls apart with a constant speed. The response force may be easily integrated over distance, giving the work required to stretch the system. The maximum work W_max required to stretch the system is also bigger if oscillations were introduced before. Again, the difference is largest for glutenins, this time followed by maize proteins. Sometimes the force drops to near zero after reaching the maximum (meaning that the system was ruptured) and the work stops rising (as seen for the dark blue curves on Fig 3). There are other cases where a substantial mechanical work is still required for elongating the system after reaching F_max (cyan curve), but the total mechanical work W_max is correlated with the height of the F_max peak. The proteins apply some effective pressure on the walls, however it is much smaller than the mechanical stress resulting from the stretching (the atmospheric pressure would exert 0.01 nN of force on the walls).

Even if entanglements cause sharp peaks in the force curves, the resistance against deformation is also controlled by the number of inter-chain contacts, n_inter, that have to be ruptured during elongation. The biggest increase in n_inter is observed for maize proteins. They have a lower number of contacts in total (as indicated by their coordination number z shown in Table 1), but their residue composition allows them to make more contacts due to the mixing effect of the oscillating deformations.

The oscillatory mixing effect is general and it operates in each of the systems studied. It results in a stronger network of proteins that are connected together more tightly, which results in the lowering of RMSF. Less mobile proteins may be more sturdy, but, as a trade-off, the maximum force may occur for a lower inter-wall distance s [51].

The molecular explanation of gluten elasticity

According to the “loops and trains” picture, the response of gluten to deformation consists of two stages: at the first stage the chains elongate and the cavities between them shrink, which requires smaller force than disentangling the chains and moving them alongside each other that takes place in the second stage [14, 52]. We tested these predictions by measuring how 6 different quantities change during the final extension of the simulation box (see Fig 4). The snapshots shown in this figure indicate that the “trains” do form, and that shearing between the chains drawn in red and blue may lead to a substantial force by disrupting and reforming the contacts between the chains. Indeed the number of inter-chain contacts decreases in the beginning of the elongation, but then it increases again as “trains” are formed. The chains become more rod-like and less globular, as indicated by the distortion parameter w [53] (averaged over all chains) that is equal to 0 for an ideal sphere, and 1 for an ideal rod. It is defined as $$ \frac{\Delta R}{\overline{R}}$$, where $$ \overline{R}=\frac{1}{2}(R_1+R_3)$$ and $$ \Delta R=R_2-\overline{R}$$; R₁ is the smallest radius of inertia, R₃—the largest and R₂ the intermediate.

Fig 4

The time dependence of 6 different quantities that are measured for gluten (as defined in section 1 of the S1 Text) during the elongation in the last stage of the simulation.

The solid blue (orange) lines correspond to the quantities listed along the left (right) y-axis. The lines are smoothed, solid (hollow) dots represent raw data from simulation. The quantities analyzed are: the number of inter-chain contacts n_inter and the average distortion parameter w, force F and work W, the maximum volume of a cavity $$ {\mathrm{V}}_{\text{max}}^{\mathrm{c}}$$ and the number of cavities n^c. The 4 snapshots (similar to those on Fig 1, but in the ortographic projection) correspond to the displacements s written above them and marked by the arrows. The same coordinate system (defined in the top left) is used for all snapshots. The biggest cavity is outlined by a pink mesh.

Formation of a “train” (two parallel chain fragments connected by hydrogen bonds) does not coincide with the maximum of the force, because F_max occurs near the beginning of the elongation and is associated with the biggest decrease in the number of contacts. Nevertheless, the force does not drop to zero and the work required to stretch the system steadily increases, because the chains from the opposite ends of the box are still connected and cause further resistance against deformation.

In order to quantify the number and size of cavities that form in gluten, we used the Spaceball algorithm [15, 16], which uses a probe with a user-defined radius to determine a three-dimensional contour of the system. The volume of the cavities is computed by filling a grid with those probes and rotating the grid with respect to the system. Our system is quite large, so we used grid with 0.1 nm lattice constant, a probe with radius 0.38 nm and one rotation of the grid. The snapshots generated show that the biggest detected cavity is dispersed over the whole simulation box, but during the elongation it becomes associated with the set of proteins attached to the Z = 0 wall (for the trajectory shown). As they are more and more elongated, the size of the cavity decreases. Larger cavities become divided into many smaller cavities, which is reflected in the increase in n_c.

Fig 4 illustrates just one simulation (the same properties for the same system, but with a different random seed, are shown in S6 Fig), but the volume of the biggest cavity indeed decreases after elongation in most cases (see S4 Fig). The chains become more rod-like not only during the final stretching, but also after the oscillations (see S2 Fig): deformations cause the proteins to be briefly stretched, but after the deformation ends, the proteins do not return to the exactly the same shape as before, but stay elongated in one direction (which may be called a memory effect).

Determining the dynamic shear modulus

One of the few fundamental (independent of the size and shape of the sample) rheological properties of gluten is the dynamic shear modulus G* = G′ + iG′′, measured by deforming the gluten sample with sinusoidal deformation (shear strain for G* and normal strain for E*) [41]. The real part (G′, elastic modulus) represents the in-phase part of stress response, the imaginary part (G′′, loss modulus) the out-of-phase response (phase difference δ is defined as $$ tan\delta =\frac{G^{\prime \prime }}{G^{\prime }}$$). This is a non-destructive method, because the amplitude of oscillations is small (usually up to a few percent [2]). The total force F(t) that the system under a shear deformation exerts on the walls in the X direction (see S1 Fig), divided by the surface area of a wall S, defines shear stress ϕ(t) = F_X(t)/S (the usual symbol for shear stress, τ, is a time unit in our model, so we use ϕ instead). A strain defined as γ = γ₀ cos(ωt) is expected to induce a stress in a similar form, shifted in phase by δ: ϕ(t) = ϕ₀ cos(ωt + δ). Dynamic shear modulus G* = G′ + iG′′ is then defined as [54]:

In Table 2 and Fig 5 we compare the results with experiment. Table 2 shows that gliadins are responsible for the viscous properties of gluten (big δ), and glutenins for the elastic part (small δ). Despite largely different values of the results, this trend is still visible in our simulations for the longest oscillation period (70 μs): tan δ_sim is biggest for gliadins and smallest for glutenins.

Fig 5

Dynamic shear modulus (top panel) and tan δ (bottom panel) as a function of oscillation frequency f.

Results are from our simulations (on the right) and from experiments (on the left). Purple points are from [57], red and light green from [56]. On the top panel, orange and dark green points are from [58]. All curves fitted to the data are of the form a ⋅ f^b + c, where a, b and c are the fit parameters.

Table 2

Shear modulus of dried gluten and its fractions, found by experiment with oscillation frequency f = 1 Hz [2], and by our simulations (marked by subscript_sim) with f ≈ 14 kHz (corresponding to the longest oscillation period we used, 70 μs).

System	G′ [kPa]	G′′ [kPa]	tan δ	$$ [MPa]	$$ [MPa]	tan δsim
Gluten	2.5	1.3	0.5	6.0 ± 1.6	6.1 ± 1.7	1.0 ± 0.1
Gliadins	0.5	0.6	1.2	5.6 ± 1.7	8.0 ± 2.3	1.4 ± 0.1
Glutenins	40	9	0.2	8.3 ± 2.2	7.9 ± 2.1	0.95 ± 0.1

We use frequencies that are 4 orders of magnitude larger than in experiment, but some extrapolations can be deduced from the experimental frequency spectra of gluten’s G*. We used the results from two experimental papers [56, 57] that use hydrated gluten samples. Both predict that tan δ should be rising with frequency—see the bottom panel of Fig 5, where tan δ_sim is between the extrapolations predicted by [56] (green curve) and [57] (purple curve). Such large discrepancies are the effects of different gluten compositions, that greatly influence the results (we use the Aubaine dataset from [57]).

Growing δ means that G′ and G′′ converge, as seen in the top panel of Fig 5. Both G′ and G′′ are rising with f, however this rise is insufficient to explain why the values from simulation are 3 orders of magnitude larger. Either there is a significant non-linearity for higher frequencies and the extrapolation no longer holds, or this is a limitation of our model. We use implicit solvent which cannot account for differences in gluten’s hydration that may greatly affect the results (even “dry” gluten is slightly hydrated) [56].

On the other hand, the results of the simulation are in good agreement with the experiments done on gluten-based bioplastics (as shown on the top panel of Fig 5), where gluten proteins are mixed with glycerol and melted at over 380 K [58]. Even for bioplastics, the tensile strength (the stress needed to break the deformed system) is less than 1 MPa [58], while the F_max divided by the wall area (which is around 100 nm²) can be more than 10 times higher (see Fig 3).