Escherichia coli Peptidoglycan Structure and Mechanics as Predicted by Atomic-Scale Simulations

Bacteria face the challenging requirement to maintain their shape and avoid rupture due to the high internal turgor pressure, but simultaneously permit the import and export of nutrients, chemical signals, and virulence factors. The bacterial cell wall, a mesh-like structure composed of cross-linked strands of peptidoglycan, fulfills both needs by being semi-rigid, yet sufficiently porous to allow diffusion through it. How the mechanical properties of the cell wall are determined by the molecular features and the spatial arrangement of the relatively thin strands in the larger cellular-scale structure is not known. To examine this issue, we have developed and simulated atomic-scale models of Escherichia coli cell walls in a disordered circumferential arrangement. The cell-wall models are found to possess an anisotropic elasticity, as known experimentally, arising from the orthogonal orientation of the glycan strands and of the peptide cross-links. Other features such as thickness, pore size, and disorder are also found to generally agree with experiments, further supporting the disordered circumferential model of peptidoglycan. The validated constructs illustrate how mesoscopic structure and behavior emerge naturally from the underlying atomic-scale properties and, furthermore, demonstrate the ability of all-atom simulations to reproduce a range of macroscopic observables for extended polymer meshes.

Derivation of Young's moduli E g and E p Under the assumption of linear elasticity, the constitutive equations for a material are known as the generalized form of Hooke's law, namely σ ij = C ijkm km (1) where σ ij is the stress tensor, km the strain tensor, and C ijkm is the fourth-order stiffness tensor.
Although C ijkm contains 81 constants, due to various symmetries not all of them are fully independent.
Two additional properties of peptidoglycan further simplify the constitutive equations. First, because peptidoglycan is very thin compared to the size of a bacterium, plane stress conditions are assumed, i.e., σ zx = σ yz = σ zz = 0. Second, peptidoglycan is orthotropic, i.e., it possesses two orthogonal planes of elastic symmetry [1]. Under these assumptions, Eq. 1, specifically its inverse, takes the form where the constants C ij are now expressed in terms of traditional physical quantities, the Young's moduli E x,y , shear modulus G xy and Poisson's ratios ν xy and ν yx . Expanding this equation gives for the strains in the x and y directions By solving the equations for E x and E y and assigning the x axis to the peptide cross-links and the y axis to the glycan strands, the individual elasticities are determined to be which is precisely Eqs. 2 and 3 in the main text.

Relationships between stress and pressure
In an MD simulation, the pressures in the x, y, and z directions within the periodic simulation box can be measured and, to some degree, also controlled. In a pure medium, by virtue of Newton's first law, the pressure and the stress in a given direction are opposite one another, i.e., P i = −σ i [2]. However, for simulations of the cell wall, the simulated volume is composed of both the hydrated peptidoglycan layer and a layer of water above and below. Thus, in these simulations, the pressures over the entire box, which are the quantities typically reported, can be decomposed into averages over the separate regions, the peptidoglycan layer (PG, its region given by its thickness t) and the water layers (wat., given by L z − t), giving where i = x, y, or z and the definition of the z-dependent pressures follows the formalism developed in the context of membrane surface tension [3,4]. Under hydrostatic conditions, P wat.
x = P wat. y = P wat. z = P wat. . Therefore, where σ i denotes the stress resultant (averaged over the cell-wall patch) in the i direction.
For a porous material such as the cell wall, the body force due to the water pressure can be eliminated through the definition of an effective stress, σ i = σ i + P wat. , for which the constitutive relations still hold. By substituting σ i in Eq. 8 one obtains which then can be solved for σ i , giving Finally, by solving Eq. 10 for i = z, P wat. can be expressed in terms of σ z , giving for the effective stresses in the x and y directions where P z is held fixed at 1 atm in all simulations, and, under plane-stress conditions, σ z = 0. Eq. 11 is identical to Eq. 4 derived for simulations of microtubules in Wells and Aksimentiev [5].
The thickness t in Eq. 11 is that of the stress-bearing portion of the peptidoglycan layer and is generally not identical to that based on the mass density. To measure t in each simulation, the lateral pressure profiles in 1-Å slices along the z axis were first calculated. The resulting profiles, which are further smoothed by averaging over of a sliding 5-Å window, display a sharp increase in the lateral pressure in the direction of the applied strain. The stress-bearing region was taken to be that for which the lateral pressure is more than 10% of the maximum and the corresponding width was used as t (see Fig. S1). Due to large fluctuations in the pressure and the limited number of frames available for post-processing, we assume the magnitude of the pressure in each individual window is insufficiently converged for direct calculation of the stress [4]. Thus, we instead use Eq. 11 along with the pressure over the whole box during the original simulation, which is measured every time step rather than every frame.

Validity of plane-stress approximation
While it is common to assume a condition of plane stress in many applications, proof of its validity is often lacking [6]. To explicitly examine this assumption, and to determine an estimate for σ z , pressure profiles were again used. Although noise is large, there is a constant non-zero stress in the z direction. This stress is independent of the applied strains σ g and σ p as shown in Fig. S2, and is present even for a completely unstrained patch. Rather than arising from specific mechanical properties of peptidoglycan, however, the observed σ z is due to an entropic pressure, i.e., the tendency of unlinked peptides to expand away from the central plane, akin to the pressure generated by an ideal gas.
The constant, but non-zero, σ z will manifest itself in two places in the previous derivations, one for the calculation of the stresses σ x,y from the measured pressures and one for the calculation of the elasticities from Hooke's Law. In the former case, the effect is simply additive (see Eq. 11). In the latter case, one must first consider a fully three-dimensional orthotropic material, for which the stiffness tensor becomes sixth order. Ignoring again the shear stresses, Eq. 1 becomes which, when multiplied out, gives for x with y and z following similarly. Thus, the term due to σ z can be treated simply as a perturbation on x and y in the equations derived for the plane-stress state. The elasticities in Eqs. 5 and 6 then become where the z dependence is accounted for by We recall from Eq. 11 that the determined strains σ g and σ p are now also dependent on σ z , making the final expressions for the elasticities where σ 0 g and σ 0 p denote the values under plane-stress conditions.
The elasticities are determined from the simulations by fixing one of the strains g or p at zero and fitting σ 0 g and σ 0 p as a function of the remaining strain. Re-expressing Eqs. 18 and 19 in this form gives Thus, although σ z may be non-zero in our simulations, because it is independent of the applied lateral strain (see Fig. S2), it will not affect the slope of the σ 0 -fit. As the elasticity is derived from this slope, the plane-stress approximation can safely be used here.
As another way to check the validity of the plane-stress approximation, we also consider the thickness of the cell-wall patch. In order to maintain an effective zero-stress state in the z direction (discounting the entropic pressure shown in Fig. S2), the thickness should respond freely to changes in lateral dimensions such that the volume remains constant [2]. More precisely, if the dimensions of the patch are given by L x , L y , and L z and changes in each length by dx, dy, and dz, then where the strain x is given by dx/L x,0 (similarly for y and z). For the volume to remain constant, leading to the equations and, thus, where i = x or y and t is the thickness of the patch.
By extrapolating the thickness as a function of the strain , the predicted thickness under planestress conditions, L z,0 in Eq. 25, can be found. For each of the five cell-wall patches examined, we measured the thickness in simulation as a function of g and p , weighting each by (1 + g ) and (1 + p ), respectively. Should plane-stress conditions hold, these weighted values will be independent of . On the other hand, if σ z affects the thickness in a strain-dependent fashion, extrapolation of a linear fit to = 0 will provide a different value of L z,0 than the average of the weighted thicknesses from all simulations. Table S2 compares these two metrics, i.e., the average and extrapolated values of L z,0 .
For all patches except avg26, the difference between the two metrics is less than 10% and, furthermore, is not consistently positive or negative. Therefore, because the thickness in almost all cases responds freely to the imposed strain in the plane, we again conclude that plane-stress conditions are a valid approximation for the cell-wall patches simulated.

Stress-strain relationships
To calculate Young's modulus, first the slope of the stress-strain relationship was determined based on approximately 25 simulations for each constructed patch, divided between those probing E g ( p = 0) and those probing E p ( g = 0). When this slope is weighted by (1 − ν pg ν gp ), the elasticity is recovered (see Eqs. 4 and 5 in the main text). The average stress in each simulation is plotted as a function of strain in Fig. S3. For avg8, approximately 170 ns of simulation in total were used to determine E p and E g due to the slow convergence of the average pressure. For the other four systems, approximately 50 ns each was sufficient.

Calculation of pore size
In order to calculate the largest effective pore size in a given patch of peptidoglycan, first a grid of points spaced 0.3 nm apart was overlaid on the plane of the layer. At each point, the maximum-radius sphere that could be inscribed without contacting either glycan strands or cross-linked peptides was determined (see Fig. S4). Because it is assumed they are flexible and, thus, would not impede diffusion through the layer, peptides that were not cross-linked were ignored. The global maximum radius was then taken as the maximum pore size for a given structure. Although some pores are clearly irregular, experimental estimates of pore size were also calculated under the assumption that the molecules diffusing through were spherical, at least in one of two presented models [7]. Numbers reported in Table 2 in the main text are time-averaged over the course of all simulations at the stated applied strain.

Strain-dependent strand addition
To examine how strain may affect the placement of new strands in a growing cell wall, a new system was constructed. Beginning with the avg17 peptidoglycan patch under an applied strain of p = 0.2, one of the strands was deleted. The system was then equilibrated under tension for 10 ns, permitting the remaining strands to retract from the newly formed gap. A completely extended strand of the same length and composition as the deleted one was then added back into this gap and its glycans were held in place while the peptides were free to move for a 3 ns simulation, i.e., the same procedure used for initial construction of the entire patch. Peptides side chains that came near each other during this simulation were cross-linked, maintaining the 50% ratio of linked to unlinked as before. Finally, the resulting new patch was simulated both under constant strain ( p = 0.2, g = 0) for 3 ns and also under no strain for 6 ns.
As seen in Fig. S5, some cross-links formed in new locations, resulting in a different patterning of the strands. One immediate observation is the formation of a larger pore (3.6 nm radius vs. 2.9 nm originally), due to the elimination of a connection between the end of the re-added strand and one to its left in the figure. However, the surrounding strands do not change dramatically. When allowed to relax under zero applied tension, the new patch's equilibrium dimensions changed slightly, becoming 19.4±0.3 × 32.4±0.3 nm 2 , as compared to 18.7±0.2 × 33.4±0.25 nm 2 originally. Thus, while the equilibrium area is identical (within less than 1%), the dimensions have changed by +3.7% and -3% in the peptide and glycan directions, respectively. This change is one consequence of the decrease in the diversity of strand-strand cross-links; whereas the original strand was connected to three others, the re-added one is only connected to two. Although not explicitly tested, a likely consequence of this change in cross-linking is a slightly lower E p .