1 Introduction

Gluten is very important in food industry because of its unique viscoelastic properties [1, 2] that are especially involved in the context of making bread and pastries. It is mostly an ensemble of proteins that are intrinsically disordered. We have recently proposed a protocol to study viscoelasticity of such ensembles through molecular dynamics simulations [3]. The protocol uses our novel coarse-grained model [4] that has been designed primarily for the intrinsically disordered proteins. Studies of the system’s shearing or extension require introducing a containing box. The simulation box we used had periodic boundary conditions in the X and Y directions and solid walls in the Z direction. Changing the location of the solid walls enabled studies of the gluten elasticity. Two different types of walls were simulated: an FCC-lattice of beads or a ,,sticky” wall that made an adaptive contact when a protein residue appeared sufficiently close to the wall. Here we present how these two approaches change the process of adhesion to the walls, and how the formation of the contacts with the wall change the conformations of the proteins. We also study the coordination number in the system, and how it affects its rheological properties according to recent theoretical predictions [5].

2 Methods

Our coarse-grained dynamic structure-based model is described in detail in Ref. [4]. It represents residues as harmonically connected beads. The chain stiffness is ensured by the bond and dihedral angle statistical potentials derived from a random coil library [6]. The excluded volume effect is ensured by the repulsive part of the Lennard-Jones potential with the minimum set at 5 Å. The non-local interactions include Debye-Hueckel electrostatics with distance-dependent permittivity [7] and Lennard-Jones attractive contacts, that are dynamically created and destroyed based on the local structure of the protein chain [4]. Each of the three possible types of contact (between chain backbones, sidechains and the mixed backbone-sidechain contact) is associated with different criteria involving beads’ distance, their coordination number and their geometry in the chain. When those criteria are met, the attractive Lennard-Jones potential is quasi-adiabatically turned on. When a contact breaks, the attraction is turned off in the same manner. Thus, a one-bead-per-residue model can capture structural details involving backbone and sidechain positions [4]. It is especially important for gluten, made from intrinsically disordered proteins that can change their structure dynamically. Because they do not have one unique conformation, we start the simulations from self-avoiding walk conformations.

The simulation protocol and the general properties of gluten are discussed in Ref. [3]. Here we summarize only the most important information about the simulated systems.

Gluten proteins can be divided into gliadins and glutenins, based on their solubility in alcohol solutions [2]. Glutenins form an elastic network responsible for gluten’s resistance to stretching, while more viscous gliadins act as a plasticizer [8].

In our simulations, gluten is represented by 6 different types of chains (3 for gliadins and 3 for glutenins), mixed in proportions reflecting the experiment [2]. The average number of residues in gliadins is 290. Glutenins are longer, with an average length of 445 residues [3]. The total number of residues in the gluten system studied is 4271. For computing the shear modulus, we made additional simulations of only gliadins (4322 residues) and only glutenins (3985 residues). We have also simulated smaller systems made from corn and rice proteins (3504 and 3503 residues total, respectively). The corn system is made from 3 different types of \(\alpha\)-zein chains and 1 type of both \(\gamma\)-zein and glutelin-2 chains. They are representative for the protein balance in maize seed storage proteins [9]. The rice system is constituted by 8 different types of glutelin chains, 4 types of prolamin chains and 1 type of globulin chains [10]. Those proteins are significantly shorter than glutenins (average length of 221 residues for corn and 208 residues for rice proteins).

We compare two different mechanisms of gluten proteins’ adhesion to the walls: FCC walls and adaptive walls. In the first case, each wall is made from two layers of beads arranged as a (111) triangular face of the FCC lattice. The lattice constant of the triangular lattice was taken to be 3 Å. For the adaptive walls, an interaction center appears only if a distance between the wall and a protein residue crosses the threshold of 5 Å. In both cases the protein-wall interaction is described by the Lennard-Jones potential with minimum at 5 Å  and the well-depth of 4 \(\epsilon\), where \(\epsilon \approx 110\) kcal/mol [11, 12]. However, in the case of the FCC walls a protein residue can be attracted by multiple wall sites simultaneously, while in the adaptive walls one residue is attracted by only one interaction center. In both cases a residue can detach from the wall and then reattach later in a different place.

We compute the dynamic shear modulus \(G^*=G'+iG''\) by fitting a function \(F(t)=F_0\cos {(\omega t + \delta )}\) to the force exerted on the walls by the sample during periodic oscillations of the distance between the walls (the frequency of the oscillations is \(\omega\)). Then \(G^*\) is calculated according to the following equation:

$$\begin{aligned} G'=\frac{F_0/S}{A/s_0}\cos \delta \qquad G''=\frac{F_0/S}{A/s_0}\sin \delta \end{aligned}$$
(1)

where S is the area of the walls, \(s_0\) is the equilibrium distance between the walls and A is the amplitude of the oscillations. An example of the fitted F(t) curve is shown on Fig. 1.

3 Results

3.1 Adhesion

Fig. 1
figure 1

The force F acting on the walls (violet) and the distance s between the opposite walls (green) during one exemplary simulation of gluten. Big letters indicate different stages of the simulation. Horizontal dotted lines indicate the equilibrium distance \(s_0\) and the distance \(s_{max}\) corresponding to the maximum force, indicated by a vertical dotted line (Color figure online)

Full size image

The simulation is divided into different stages: in the first stage (S), the proteins are squeezed together into a cubic box, from a rarefied phase, to reach the target density. This stage is discussed in Ref. [3]. It involves walls that are purely repulsive during squeezing. The walls are perpendicular to the Z axis (and periodic boundary conditions in the X and Y directions are used during the whole simulation). Then, after some time for equilibration, the attractive interactions with the walls are turned on (stage W). After another equilibration (stage E) the box undergoes temporary periodic deformations with the period of 20 \(\mu\)s (stage O). In the case of pulling, the deformation is accomplished by moving the walls in the Z direction. In the case of shearing—along the X direction. The next stages are the final equilibration (E) that is followed by pulling the walls apart in the Z direction (stage P), resulting in tearing the system into two subsystems. All the stages are depicted in Fig. 1. When we refer to the variables computed ,,before” and ,,after” the oscillations, we mean averages taken over the equilibration stages E preceding and succeeding the oscillation stage O.

The most important notion in our model is the concept of a contact: an attractive interaction that may arise between two beads. The beads represent protein residues and may represent sites on the rigid walls. The coordination number \(n_c\) is defined as the mean number of contacts made by a protein residue. Analogically \(n_w\) denotes the mean number of wall contacts (FCC or adaptive) per residue. Fig. 2 shows how these two variables change during different stages of the simulation.

During the initial equilibration, \(n_c\) fluctuates, because the chains are reorienting themselves.Footnote 1 Turning on the wall contacts (indicated by a sharp increase in \(n_w\)) results in lowering \(n_c\) in the case of FCC walls, because each residue can make only a limited number of contacts [4]. Adaptive walls have lower \(n_w\), so the drop in \(n_c\) is smaller.

Before the onset of the periodic oscillations, the walls are made to move to the maximal amplitude of the oscillations, which is 10 Å (this moment is marked by letter A in Figs. 1 and 2). After moving the walls, \(n_w\) rises significantly in the case of shearing deformation and adaptive walls: shearing brings more residues close to the walls and more adaptive contacts are made. In the case of pulling, the box is squeezed to its minimal size. However, it is not accompanied by such a big rise in \(n_w\) - shearing makes more significant changes in the adhesion. Periodic oscillations (O) are reflected in the oscillating \(n_c\) (especially for adaptive walls in pull mode). It is interesting to note that the final pulling (P) does not lead to a large decrease in \(n_c\), because inter-chain contacts are replaced by intra-chain contacts [3].

Fig. 2
figure 2

Time dependence of the coordination number \(n_c\) (violet) and the number of wall contacts per residue \(n_w\) (green). The orange color marks the temporal segments in which various box deformations take place: adjusting the amplitude (A), periodic oscillations (O) and final pulling (P). The panels correspond to simulations with the FCC walls (left) or with the adaptive walls (right), in shearing (bottom) or pulling (top) modes. The lines were artificially smoothed out (the raw data were saved every 0.1 \(\mu\)s) (Color figure online)

Full size image

The change in \(n_w\) induced by the box deformation is summarized in Fig. 3. The number of contacts rises in each case, but it is significant only in the case of shearing deformation and adaptive walls.

Fig. 3
figure 3

The average number of the wall contacts per residue \(n_w\) before (violet) and after (orange) periodic deformations of the box (Color figure online)

Full size image
Fig. 4
figure 4

Snapshots from the simulations with the FCC (left panels) or adaptive walls (right panels), before (top panels) and after (bottom panels) periodic shearing oscillations. The big spheres represent centers of mass of the protein chains. Chains are represented in the balls-and-sticks convention. Each residue is colored according to its coordination number \(n_c\) and scaled from 2 (dark blue) to 10 (dark red). The dashed line on the bottom left panel outlines a rift where no covalent bonds are present (Color figure online)

Full size image

The box deformation induces changes in the shapes of the proteins. This shape can be quantitatively described by two parameters: the asphericity w and the radius of gyration \(R_g\). They are defined by [17]:

$$\begin{aligned} R_{g}&= |\frac{1}{n} \sum _{k=1}^n \mathbf {r}^{\;2}_{k}| \end{aligned}$$
(2)
$$\begin{aligned} w&=\frac{\Delta R}{\overline{R}} \end{aligned}$$
(3)

Here, n is the number of residues in the protein, \(\mathbf {r}_{k}\) are their position vectors with respect to the protein’s center of mass. The parameter w depends on all three main radii, \(R_{\alpha }\), associated with the eigenvalues of the tensor of inertia \(D_{\alpha }\): \(R_{\alpha }=\sqrt{D_{\alpha }/n}\). These eigenvalues, and thus the radii, now correspond to the instantaneous (as opposed just to the native) shape of the protein. \(R_1\) is the smallest radius, \(R_3\) – the largest, and \(R_2\) – the in-between one; \(\overline{R} = \frac{1}{2}(R_{1} + R_{3} )\) and \(\Delta R = R_{2}-\overline{R}\). Globular shapes correspond to w being close to 0. Elongated cigar-like shapes yield substantial positive values of w because then \(R_{2}\) is close to \(R_{3}\) and \(w \sim \frac{1}{2}(R_{2}-R_{1})\). Substantial negative values of w indicate planar shapes as then \(R_{2} \sim R_{1}\) and \(w \sim \frac{1}{2}(R_{1}-R_{3})\).

Figures 5 and 6 show these parameters before and after the deformation. Proteins become more elongated (their w and \(R_g\) increase) and their centers of mass move closer to the walls. In the case of FCC walls, this creates a rift in the middle of the simulation box, as visible in Fig. 4. In simulations with the adaptive walls, the chains in the middle of the box span both walls.

Fig. 5
figure 5

The asphericity parameter w as a function of the Z coordinate of the corresponding chain’s center of mass. Each point represents a different snapshot from a simulation. Snapshots were taken before (violet) and after (orange and yellow) periodic box deformation. Panels show simulations with FCC walls (left) or with adaptive walls (right), in shearing (bottom) or pulling (top) modes (Color figure online)

Full size image

Gluten proteins can be divided into large glutenins that create the network responsible for gluten’s elasticity, and smaller, more viscous gliadins [2]. This distinction is visible in Fig. 6, where the bottom of each panel represents gliadins with small \(R_g\). Bigger glutenins tend to stick closer to the walls.

Fig. 6
figure 6

The radius of gyration \(R_g\) as a function of the Z coordinate of the corresponding chain’s center of mass. Each point represents a different snapshot from a simulation. Snapshots were taken before (violet) and after (orange and yellow) periodic box deformation. Panels show simulations with FCC walls (left) or with adaptive walls (right), in shearing (bottom) or pulling (top) modes (Color figure online)

Full size image

3.2 Shear Modulus and Critical Strain

Gluten can be considered a cross-linked polymer network, because it can form covalent disulfide bonds [2] and entanglements [8]. Each contact can also be considered as a link: the coordination number is a measure of connectivity in the system. Networks that have more connections have higher elasticity [5], which is reflected in their dynamic shear modulus. For viscoelastic substances like gluten this modulus is a complex number, where the real part \(G'\) describes the in-phase elastic response to periodic deformations, and the imaginary part \(G''\) is the out-of-phase viscous part [18]. Therefore, \(G'\) should rise with \(n_c\). The top panel of Fig. 7 confirms it, although the data are very noisy: different \(n_c\) come from different systems (rice, corn, glutenins, gluten and gliadins) that have many other differences besides \(n_c\). For example, corn proteins have the possibility to create more disulfide bonds than rice proteins [19], they also have more hydrophobic residues [20], which may lead to their higher shear modulus [21]. The \(n_c\) values are averages taken from the equilibration stages before the oscillations.

Another property that changes with \(n_c\) is the critical strain: for higher strains the sample is torn apart. This is realized in our simulations by recording the distance \(s_{max}\) between the walls, for which the highest maximum in the force occurs during pulling apart the walls [3]. The distance in equilibrium is \(s_0\), so critical strain \(\gamma =(s_{max}-s_0)/s_0\). Examples of \(s_0\) and \(s_{max}\) are shown on Fig. 1, while the critical strain \(\gamma\) is shown on the lower panel of Fig. 7. Counterintuitively, lower \(n_c\) leads to higher \(\gamma _0\). This is explained by the slip-bond theory [5], which predicts that networks with lower connectivity can re-form the contacts more easily during breaking, and networks with higher connectivity are more susceptible to a critical rupture.

Fig. 7
figure 7

The real part of dynamic shear modulus \(G'\) and the critical strain \(\gamma _0\) as a function of coordination number \(n_c\) for the 5 studied systems (top labels). The gluten system is the same as in the previous figures

Full size image

4 Conclusions

When shearing is induced in simulations with adaptive walls, the proteins exhibit much greater increase in \(n_w\) (which correlates with the adhesion strength) than in the FCC case: for ,,rough” FCC walls each contact with a new bead is compensated by a loss of a contact with a bead that moved further away, thus \(n_w\) rises only slightly. Contrary, each residue can make only one adaptive contact with the wall: if it is broken, it can only be reformed later. The large increase in \(n_w\) for the adaptive walls is possible because the total number of wall contacts is smaller than in the FCC case. This allows the protein chains to move further away from the wall to also make contacts with it. Those ,,central” chains become more elongated and act as ,,transmitters” of the strain in the system. This leads to much better rheological properties when compared to the FCC walls (see Fig. S10 in the Supporting Information of Ref. [3]). In the FCC walls, each protein can make a large number of contacts with a wall, which leads to dense boundary layers near the walls and more dispersed and viscous layer in the middle (see Fig. S9 in the Supporting Information of Ref. [3]).

What is described above corresponds essentially to the same effect that is predicted by the slip-bond theory [5]: networks with higher \(n_c\) can be more susceptible to rupturing, because when strongly interconnected parts are moved, they are too stiff to adapt to the changes so they can no longer keep the connections. On the other hand, more loosely connected networks have the ability to re-form the contacts and quickly adapt. Our results are consistent with this picture, not only in the case of critical strain, but also in the context of adhesion to the walls: FCC walls make more contacts (higher \(n_w\)), but adaptive walls have better rheological properties [3]. Instead of forming two layers closely connected to their respective walls (like FCC walls), adaptive walls allow slipping and reforming of contacts, so the system is more homogenous and does not divide into layers. This novel application of the slip-bond theory [5] to wall adhesion during shearing could be in principle replicated in experiments by measuring distances of the protein from the stiff walls.