Molecular Dynamics of Mycolic Acid Monolayers
Mycobacterium tuberculosis, the organism responsible for TB infection in humans, is inherently resilient against host defences and anti-TB drugs. This persistence is attributed partly to the presence of lipids, such as mycolic acids (MAs), which make the cell wall impermeable. To study the conformational dynamics of MAs, we present a coarse-grained model for a representative α-MA (AMA) from Mycobacterium tuberculosis using the MARTINI force field. The model is used to simulate monolayers of different sizes; a small monolayer consisting of 220 MAs and a large monolayer consisting of 1972 MAs. The model could replicate key features of experimental monolayers such as phase changes and the collapse point. By studying the conformation of MAs in the simulated monolayers, it was 1 found that AMA did not fold into the W-conformation at large surface areas but was only folded at the head group to give a wide U-shape. On monolayer compression, the MA chains came closer together, into a narrower U-shape, and an ordered monolayer was formed before it collapsed.
Introduction<!>Methods<!>Atomistic simulations<!>Coarse grained simulations<!>Parameterisation of bonded interactions<!>Monolayer setup<!>Monolayer simulations<!>Surface tension reference values<!>Analysis details<!>Optimisation of bonded parameters<!>Alpha-MA isotherms<!>Mycolic acid conformation in monolayer compression<!>Conclusion
<p>Despite being largely a curable disease, 1.6 million people died from tuberculosis (TB) in 2017. 1 There are numerous factors that complicate the control of TB. These include co-infection with human immunodeficiency virus (HIV) as well as development of drug resistant strains. 1 In addition, long treatment regimens (at least 6 months) on a combination of drugs with severe side-effects leads to patient non-compliance. The sensitivity of some diagnostics (especially in HIV-co-infected patients) is low and the long time required for confirmed diagnosis (2-6 weeks) plus an additional 4-6 weeks for resistance determination by bacterial culture 2,3 complicates matters further.</p><p>Mycobacterium tuberculosis, the organism causing TB infection in humans, is inherently resilient against the host defences and anti-TB drugs. This is due, in part, to its lipid-rich cell wall making it impermeable to many currently employed antibiotics. Mycolic acids (MAs) are β-hydroxy acids of up to 90 carbons long and make up part of the outer layer of the cell wall. MAs are covalently attached to arabinogalactan where they are thought to form an impermeable monolayer arrangement, but they also occur as part of the outermost layer of unbound lipids as sugar-esters.</p><p>There are three main types of MAs in M. tb, as illustrated in Figure 1A. All three types share a common 24-carbon alpha chain and β-hydroxyl acid with R,R stereochemistry. The distinction between the three types of MA is in the 2 functional groups in the mero-chain (Figure 1A X, Y). Alpha-MA (MA), the most abundant MA type in M. tb, has cis-cyclopropane groups at both the proximal and distal functional groups (Figure 1A X and Y, respectively) whereas methoxy-and keto-MAs may have either cis-or -methyl trans-cyclopropane at the proximal position (Figure 1A, X) and a methoxy-methyl or keto-methyl group in the distal position, respectively (Figure 1A, Y). 4 Figure 1: The general structure of MAs found in M. tb (A), coarse grained beads, and bead types mapped onto the all-atom structure of AMA with points a-e indicated (B) and the names of coarse grained beads for AMA (C).</p><p>The precise structural and functional roles of the different MA types are only beginning to be unravelled. It has been shown that different MA types have different effects on the immune system [5][6][7][8][9][10][11][12][13][14] and that changing the MA composition alters the outcome of an infection as well as the cell wall permeability. [15][16][17][18][19][20][21] Despite these advances, it is still uncertain how MAs are arranged in the cell wall. Zuber et al. 22 established that the bilayer of the outer membrane in which MAs are found is 7-8 nm thick. Some of the longer MAs would not be able to fit into this space without folding or bending. Therefore, the authors proposed a zipper model, in which MAs would not only have to fold at each of the functional groups, but also intercalate with lipids in the outer leaflet in order to fit into this space. When the MA is folded at all three functional groups, it creates a W-shape when mapped out in two dimensions, as shown in Figure 2 (W). We have previously defined a range of folds accessible for MAs from M. tb that are based on folding at the functional groups, 23 as well as more globular conformers. 24 A monolayer of MAs on water has been suggested as a suitable model for the tightly packed MA layer bound to arabinogalactan in the cell wall. Several such monolayer experiments have been performed [25][26][27][28][29][30][31][32][33] and some include a range of different temperatures. [29][30][31][32][33] The general trend in all these experiments is for MAs to have a folded conformation at large surface areas and to unfold as the monolayer is compressed, until the mero-chain is fully extended (with the exception of keto-MA, which stays folded, even at high lateral pressures). [25][26][27][28][29][30][31][32][33] Figure 2 shows a schematic of the unfolding that is expected in a MA monolayer. Starting from a folded W-conformation, which has a large surface area with four parallel chains, the unfolding is depicted with the assumption that a bend at the head group (b) is necessary to maintain contact of the head group with the water surface in the monolayer setup. Therefore, with the head group staying in contact with the water, an extension at the proximal cyclopropane group, c, seems the most rational option on compression. This results in a shape that is folded at the head group (b) and the distal cyclopropane group (d), which will be referred to as sZ. 23 The Z-fold represents a molecule folded at two functional groups with three chains in parallel, and sZ has symmetry in its twodimensional folding representation. Although an aZ fold, with the a point at its extended chain terminus, is also a possible conformation it is not included here as it is assumed that the molecule must bend at position b for the head group to be in contact with the water interface. It is also possible for the molecule to extend at the distal cyclopropane group, staying folded at groups b and c (eZ, with its extended chain terminating with e) as indicated by the dashed arrow in Figure 2. However, this is thought to be unlikely, as the mero-chain would not have space to extend past the water interface. On further compression it is expected that the mero-chain will extend fully so that the molecule now has a U-shape, only bending at one functional group, having two parallel chains. In a monolayer, eU is the only sensible U-shape, bending at the head-group with the tail point e at the extended chain terminus.</p><p>The eU conformation only has two parallel chains, hence occupying a smaller surface area. Any further extension in the molecule would lead to a straight MA, which is fully stretched out without bending at any point. A straight conformation is not expected in a monolayer while the head group is assumed to be at the water interface.</p><p>Due to the very long carbon chains, MAs are extremely hydrophobic and hence insoluble in most solvents. Being waxy, they also do not crystallise, making it very difficult to get structural data experimentally. For this reason, here we explored MA conformations and packing effects through molecular dynamics simulations, a powerful tool to study cellular membranes. [34][35][36] The monolayer data available from experiments [25][26][27][28][29][30][31][32][33] are used for comparison, making this an ideal system to model. Previous modelling of MAs is largely on single molecules. 23,24,30,31 These models give valuable insight to the preference of different MAs for the W-conformation, with the -methyl trans-cyclopropane group in oxygenated MAs promoting tight folding into a W-conformation more than a cis-cyclopropane group. 32 MA structure has further been shown to determine conformational behaviour. 23,24 Single molecules in vacuum, polar and nonpolar solvents were seen to adopt orderly folded conformations such as the W-conformation and amorphous globular conformations. 23,24 However, single molecules are not an accurate representation of MAs in their natural cell wall environment. MAs have also been modelled in cell wall components in which MAs were arranged as a monolayer. 37 However, in these models, MAs were truncated to have merochains of the same length as the alpha-chains, and hence these models could not show how MAs would fold and pack in the cell wall environment.</p><p>Here we apply the popular MARTINI 38,39 coarse grained (CG) force field to create a model for AMA, the most abundant MA in the M. tb cell wall. The model is used to construct AMA monolayers of different sizes. The MARTINI model was originally designed for the modelling of lipids and has been used to successfully model various lipid systems. 39 It has also been shown to reproduce lipid monolayer properties well, although quantitative reproduction of experimental isotherms remains challenging. [39][40][41][42] The use of a CG model is warranted due to the slow anticipated molecular rearrangements upon monolayer compression. Working with CG resolution has the advantage of having fewer interactions to calculate, making it more feasible to simulate large systems for longer timescales. However, the explicit stereochemistry of the functional groups in the MA is lost upon coarse-graining and it remains to be seen how this affects the accuracy of the model.</p><!><p>Molecular dynamics simulations were performed using Gromacs 4.5.5 43 and images were rendered using VMD. 44</p><!><p>Atomistic simulations of 100 ns were performed on a single AMA in vacuum, water and hexane to provide reference distributions of bonds and angles for parameterisation of the coarse-grained model. The atomistic AMA structure used is shown in Figure 1B. The stereochemistry for the chiral centres, are R,R for the acid head group and S,R for both cyclopropane groups, from left to right, as in Groenewald 2019. 24 For atomistic simulations of AMA, the OPLS all atom (AA) force field 45 (with additional parameters for cyclopropane as described in Groenewald et. al 24 ) was used and each molecule was placed in the centre of a 10 × 10 × 10 nm box and minimized with a steepest descent algorithm with a maximum stepsize of 0.01 nm, a maximum number of 200 000 steps and a tolerance of 10 kJ.mol -1 .nm -1 .</p><p>The equation of motion was integrated using a leap-frog algorithm 46 with a time step of 1 fs and the neighbour list was updated every 10 fs. Van der Waals interactions were modelled by using a shift function between 0.8 and 0.9 nm and electrostatic interactions were modelled by using PME. 47,48 V-rescale temperature coupling 49 was used at 300 K with a time constant of 0.5 ps and no constraints were applied. The simulations were 100 ns long. Vacuum simulations were performed with a NVT ensemble.</p><p>TIP4P 50 water was used to fill the simulation box for simulations in aqueous solution. Equilibration was performed with position restraints by applying a force of 1000 kJ.mol -1 .nm -2 in the x-, y-and z-directions on all carbon and oxygen atoms of the MA. The system was equilibrated by performing 100 ps NVT, followed by 50 ps NPT using Berendsen pressure coupling 51 to scale the box in an efficient way at the beginning of the simulations, and lastly 100 ps NPT using Parinello-Rahman pressure coupling 52,53 in order to yield the correct ensemble. In both NPT-ensembles isotropic pressure coupling at a pressure of 1 bar was used with a 1 ps time constant. Production simulations for MAs in water also used the Parinello-Rahman setup.</p><p>A hexane solvent box was built and equilibrated before addition to MA simulation boxes. A 3.6 × 3.6 × 3.6 nm box was filled with 216 hexane molecules, minimised (as before) and equilibrated at 300 K using a NVT ensemble for 5 ns.</p><p>The density of the hexane box was 584.201 kg.m -3 , which is approximately 11 % lower than the experimental density for hexane. 54 The simulation box containing MA was filled with hexane using the equilibrated hexane box. Equilibration was performed with position restraints by applying a force of 1000 kJ.mol -1 .nm -2 in the x-, y-and z-directions on all carbon and oxygen atoms of the MA. The system was equilibrated by first increasing the timestep from 10 -7 ps to 10 -3 ps in 10-fold increments by running respective simulations of 100 000 steps at 0 K.</p><p>Secondly, using a 10 -3 ps timestep, the temperature was increased in 50 K increments from 0 K to 250 K with consecutive 100 ps simulations and a 1 ns simulation at 300 K. Simulations were run using Berendsen pressure coupling at a pressure of 1 bar and a time constant of 1 ps. Production simulations in hexane were run using a timestep of 2 fs and with all bonds constrained.</p><!><p>The coarse-grained AMA model A mapping of 3 to 4 heavy atoms to one CG bead was used (Figure 1B). All beads were standard MARTINI 38 C1 (grey), P2 (orange) and P3 (red) types (Figure 1, B and C). Cyclopropane beads, M1A and M1B (green) are also C1 beads.</p><p>A single coarse-grained AMA was obtained from the atomistic model by using the g_fg2cg program that is part of a modified version of Gromacs, by specifying the mapping of atoms to coarse grained beads and then performing a minimisation. The single AMA was minimised with a steepest descent algorithm with an initial stepsize of 0.01 nm, a maximum number of 50 000 steps and a tolerance of 10 kJ.mol -1 .nm -1 . This conversion from an atomistic structure followed by a minimisation was performed for the starting structures of single MA simulations used in parameterisation, as well as for the reference W, sZ, eZ, eU and straight conformers.</p><!><p>Only bonded interactions were adjusted to resemble reference bond and angle distributions from mapped atomistic simulations in vacuum, water and hexane as accurately as possible. Dihedral angles were not explicitly added to the model. Atomistic simulations were performed as outlined above.</p><p>For the CG simulations on single AMAs for model parameterisation, 100 ns simulations with the standard parameters for the MARTINI 2.0 force field were used. Lennard-Jones interactions were shifted to zero between 0.9 and 1.2 nm and the Coulomb potential was shifted to zero between 0 and 1.2 nm. The relative dielectric constant was 15 and a stepsize of 20 fs was used in all CG simulations. CG simulations in vacuum were run using a leap-frog stochastic dynamics integrator. 55 Butane was used as CG solvent comparable to the hexane atomistic solvent. After water or butane solvents were added to the simulation box containing one AMA, the system was minimised using a steepest descent algorithm with an initial stepsize of 0.01 nm, a maximum number of 50 000 steps and a tolerance of 10 kJ.mol -1 .nm -1 . For temperature and pressure, the MA and solvent were coupled to Berendsen pressure-and temperature baths. The simulation temperature was 300 K with a time constant of 2 ps and isotropic pressure coupling with a 3 ps time constant, a compressibility of 3 × 10 -4 bar -1</p><p>and reference pressure of 1 bar was used for simulations with solvent. The complete set of bonded parameters is given in Supplementary Table S1.</p><!><p>The monolayer setup consisted of a central slab of water with an MA monolayer on either side, separated from periodic images in the z-direction by vacuum (shown on the left in Figure 6). Monolayers were constructed using the membrane-generating program Insane. 56 MAs were inserted into the monolayer in an approximate W-fold. Two monolayer setups were used.</p><!><p>Initially, the monolayers were minimised and then equilibrated at constant volume and temperature. Minimisations were performed with a steepest descent algorithm with an initial stepsize of 0.01 nm, a maximum number of 50 000 steps and a tolerance of 10 kJ.mol -1 .nm -1 . The small monolayer was equilibrated for 1 μs and the large monolayer for 2 μs before being compressed. The standard parameters for the MARTINI 2.0 force field were used for these simulations as described above. The temperature was set to 300 K using a Berendsen thermostat with a coupling constant of 2 ps.</p><p>The monolayers were compressed and allowed to equilibrate in a stepwise manner by alternating compression simulations with constant volume and temperature simulations, to allow the MAs to rearrange. Simulations at constant volume and temperature were run with the parameters as set out for the equilibration simulations.</p><p>Compressions were performed by setting the surface tension to -100 bar.nm (-50 bar.nm per surface) in the x-and y-directions and a compressibility of 5x10 -5 bar -1 . The compressibility in the z-direction was set to zero to keep the box size fixed in the direction normal to the monolayer. Berendsen pressure coupling was used with a coupling constant of 2 ps. The simulations under compressive stress were run until a surface reduction of ~10% was reached.</p><p>For the small monolayer, subsequent simulations at constant volume and temperature were initially 1 μs long. From A MA = 0.79 nm 2 and smaller, the constant volume and temperature simulations were increased to be 5 μs long, allowing more time for the molecules to rearrange, and for fluctuations in the surface tension to decrease in the monolayer. The initial surface tension of -50 bar.nm per surface correlates to a lateral pressure of ~2 bar. The monolayer became harder to compress at lower surface areas (even up to 10 μs simulation would not yield the 10 % decrease in surface area) and therefore the surface tension was progressively increased to -2500 bar.nm and beyond, as required, to yield the 10 % lower surface area within 1 μs simulation. The surface tension used for each compression, and the time to reach the lower surface area are listed in Supplementary Table S2. After reaching A MA = 0.23 nm 2 a sufficient decrease in surface area could no longer be achieved during compression, even with surface tension up to -4000 bar.nm (-2000 bar.nm per surface).</p><p>The large monolayer was simulated in the same stepwise approach. Constant volume and temperature simulations were performed for 2 μs before the next compression. Once a reduction in surface area of ~10% was not obtained in a 1 μs simulation, the surface tension was adjusted to -2500 bar.nm, as was used in the small monolayer. This was sufficient for the remainder of the compressions performed. The surface tension used for each compression is listed in Supplementary Table S2.</p><!><p>To obtain reference values of surface tension for water-vacuum, lipid chainvacuum, and water-vacuum, additional systems were simulated by modifying the small monolayer setup. To obtain the water-vacuum surface tension, all MAs were removed, and a 50 ns simulation was run with NVT. For the lipid chainvacuum surface tension, the lipid head groups (P2 and P3 beads) were removed and the lipid molecule linked through the C1 beads. The simulation was run at NVT. Beads at the chain termini as well as at each of the three folding points were position restrained with a force of 1000 kJ.mol -1 .nm -2 in the z-direction, so that the monolayers do not merge. For the lipid chain-water surface tension, the box used in the lipid chain-vacuum simulation was shortened to 15 nm in the zdirection and filled with water. The simulation was run for 50 ns using semiisotropic coupling with a pressure of 1 bar in the z-direction and zero compressibility for the x and y direction to keep the surface area fixed.</p><!><p>Surface tension, principal component analysis, average structures, order parameters and density profiles were calculated from the last 500 ns of the simulation at constant volume and temperature at the relevant A MA . All MAs (from both monolayers, either side of the water slab) were included in the conformational analysis.</p><p>The last 500 ns of simulation at each A MA was divided into five 100 ns blocks (excluding the last frame, so that the blocks have an equal number of frames) for plots of surface tension, surface pressure and compressibility. Block averages were obtained, from which the average is plotted, and the error bars represent the standard deviation.</p><p>The relationship between the pressure and surface tension is given by γ m = (P N -P L ) × L Z /2 where P N and P L are the normal and lateral pressures of the box with P L = (P xx + P yy )/2 and L z is the box size in the z-direction.</p><p>The compressibility modulus for the systems simulated was calculated as</p><p>All principal component plots in this work utilise the same principal components. To create the principal component space onto which to project MA conformations from simulations, in relation to reference structures, the relevant reference structures (W, sZ, eZ, eU, and straight) were converted to coarse grained molecules from carefully selected atomistic models, minimised as described before and relative energies of conformers calculated. Using these five reference structures, the eigenvectors were created and analysed by fitting to the straight conformation. Reference conformers were plotted onto the resulting 2dimensional projection of the first and second principal components, to display the most variation. Principal component analysis for the large monolayer was performed on the AMAs from every 100 th frame (every 2 ns) of the last 500 ns of the NVT simulations, after correcting the periodic boundary conditions and extracting the individual MAs. The individual MAs were also used as input for the average structure at the relevant surface area.</p><p>The uniaxial order parameter, S, was calculated for all the bonds between the coarse-grained beads with the formula: S = 0.5*(3*<cos²θ> -1), where θ is the angle between each bond and the monolayer normal and averaging is done for all the MAs over the last 500 ns of the simulation.</p><!><p>Resemblance of CG distributions to vacuum and hexane atomistic distributions were prioritised over those in water, since they more closely approximate the environment in the monolayer simulation where AMA will be at the watervacuum interface and each molecule will be adjacent to other identical, non-polar molecules. Bond force constants for CG AMA were higher than usual for MARTINI, since distributions were narrow -especially around the cyclopropane and two polar head group beads. In general, the angle force constants were also higher than usual for MARTINI beads.</p><p>Overall, bond and angle distributions correlated well with the mapped atomistic distributions, as shown in Figure 3, Supplementary Figure S1 and S2.</p><p>In particular, the distributions involving the cyclopropane beads (M1A, M1B) and polar head group (COH) match the AA distributions closely (Figure 3). This is important for a realistic description of the folding behaviour of MAs in compressed monolayers.</p><p>The angle OOH-COH-C1D did not have a very good resemblance to the reference distributions (Figure 3E). The same is apparent from dihedrals involving the COH functionality (Figure 4) This compromise had to be made to gain better resemblance of the other bonds and angles including beads COH, OOH and C1D. The bonded interactions involving the carbon chain that would determine AMA conformation were favoured over the peripheral OOH-COH-C1D angle because this angle may be less influential on the AMA conformation, with the P1 and P2 beads likely to be in contact with the water beads in the monolayer.</p><p>The two cyclopropane beads, M1A and M1B, were treated as identical in this model, with identical adjacent beads (including mapping of the beads) and identical bond and angle parameters. The CG model reproduces the atomistic angle distributions very well in hexane. However, the angles that centre on the M1A and M1B beads are represented less accurately in vacuum and in water (Figure 3, C and D) where it can be seen that smaller angles are more likely for the distal cyclopropane (M1A) than the proximal (M1B), and this difference is missed by the CG distribution. Whether this shortcoming alters the conformational behaviour of the individual molecule, and hence of the monolayer properties, is still to be determined. Smaller angles are excluded from the CG distributions and hence folding at the cyclopropane beads may be underrepresented -more so for the distal group, M1A.</p><!><p>Two monolayer systems were simulated: a small system containing 220 MAs (Figure 6) and a large system containing 1972 MAs (Figure 7).  The surface tension-area isotherms are shown in Figure 8A. For the small monolayer at a large surface area, a monotonic decrease in surface tension is seen. This is expected behaviour due to better packing of the MAs on the water/vacuum interface. A small plateau region, likely caused by the coexistence of LE and LC phases, follows around AMA = 0.6 nm 2 . A further drop in surface tension is seen below the plateau when the monolayer is in the LC phase. A shift to negative surface tension is observed at AMA < 0.46 nm 2 due to headgroups being expelled from the interface. Eventually, full collapse of the monolayer is seen. The large monolayer shows the same trends, albeit with a more pronounced LE and LC coexistence plateau region, and a release of compressive stress at small surface areas due to monolayer buckling.</p><p>During compression of the small monolayer, at around 0.79 nm 2 , some head groups of AMA molecules started to move away from the water, forming protrusions into the membrane (Supplementary Figure S3). From this point on to smaller A MA , it is seen that the longer the simulation is carried out at fixed surface area, the higher the surface tension obtained. This correlates to the MA molecules forming protrusions, breaking away from the water surface and moving toward the vacuum interface causing the molecules at the water interface to be more spread out again, obtaining a higher surface tension.</p><p>MAs leaving the water interface early on result in a difference between A MA plotted on isotherms and the true A MA . This difference becomes larger at smaller A MA as more MAs are forced away from water. For the small monolayer, it is estimated from density distributions of the head groups (Supplementary Figure S4) that the error is 0.0075 nm 2 per MA at A MA = 0.68 nm 2 , increasing to 0.1 nm 2 per MA at A MA = 0.23 nm 2 . adds stability. Alternatively, it may imply that a bilayer is preferred, as in the cell wall where free lipids effectively form a bilayer with the MA monolayer and may help to stabilise the MA monolayer. This observation could also be an artefact of our model or of limited sampling, which we are unable to determine.</p><p>No collapse plateau is seen for the small monolayer system at 300 K, likely because the monolayer is too small to bend significantly, and hence molecules are not displaced out of the monolayer plane (Figure 6). A drop to negative surface tension is seen at A MA = 0.34 nm 2 , indicating that the monolayer is unstable. The subsequent increase in surface tension would be due to a new monolayer forming at the water surface after numerous MAs have been pushed to the vacuum.</p><p>Each point on the surface tension-area plot for the large monolayer system in Figure 8A is obtained from the last 500 ns of a 2 μs simulation. MAs are also seen to leave the water interface in the large system, starting with very few molecules at 0.91 nm 2 and more at smaller A MA . The large monolayer system showed a collapse point at 0.41 nm 2 followed by a plateau. The monolayer was large enough to accumulate collective out-of-plane displacements. The monolayer starts to bend slightly at 0.71 nm 2 and has very large undulations at smaller values for A MA such as 0.30 nm 2 (Figure 7). Monolayer bending in addition to MAs moving away from the water interface results in larger differences between A MA plotted on isotherms and true A MA in the large monolayer. The approximate error in A MA due to MAs moving away from water, as estimated from density distributions of MA head groups (Supplementary Figure S5), is 0.001 nm 2 per MA at A MA = 0.91 nm 2 , increasing up to 0.14 nm 2 per MA at A MA = 0.41 nm 2 .</p><p>After the collapse point errors are estimated to be as large as 0.5 -0.8 nm 2 per MA, which is in line with numerous MAs being expelled from the monolayer at collapse.</p><p>Experimental isotherms are reported as surface pressure (Π)-area plots. To convert the surface tension (γ m )-area plots so that they are comparable to experimental data, the surface tension for the lipid chain-water, lipid chainvacuum and vacuum-water interfaces were obtained for the CG model. The surface tension obtained for lipid chain-water (43 mN.m -1 ) and lipid chainvacuum (14 mN.m -1 ) compare reasonably well with experimental values for alkanes (52 and 22 mN.m -1 , respectively) 38 and values reported for MARTINI. 31 However, the surface tension of 32 mN.m -1 obtained for the water-vacuum interface, γ aw , is significantly lower than the experimental value of 72 mN.m -1 . 38 Due to this difference, the usual conversion from surface tension to surface pressure, Π(A MA ) = γ awγ m (A MA ), cannot be applied as it stands. A simple modification in the form of Π exp (A MA ) = γ aw *γ m (A MA ) is applied to correct for the difference in surface tensions, with γ aw * a constant. 31 When relating the curve from simulation to experimental phase transition surface pressures or collapse surface pressures, it was found that the surface pressure at larger A MA is negative. This is an unrealistic scenario as it implies that the surface tension with the MA monolayer is larger than the water-vacuum surface tension, while it is usually lowered. Therefore γ aw * = 63.12 mN.m -1 was applied to the large monolayer system and γ aw * = 66.43 mN.m -1 to the small monolayer system (Figure 8B) to produce positive surface pressure values.</p><p>With this adjustment, the collapse surface pressure of 57.5 mN.m -1 for the large system correlated well with the experimentally reported value of about 60 mN.m -1 as reported by one group. 29,31 In addition, the A MA from the simulation, 0.41 nm 2 , corresponds well to the experimental range of 0.35 -0.40 nm 2 for collapse of the monolayer.</p><p>The compressibility modulus for the systems simulated, calculated as K A = A MA × δγm/δA MA , is shown in Figure 8C. The monolayers are soft at large A MA while the monolayer is in the LE phase and become stiffer as the area is reduced. Around A MA = 0.6 nm 2 , a softening of the monolayer is observed again, which corresponds to the expulsion of headgroups from the monolayer and the coexistence of LE and LC phases. This is followed by a sharp increase in the compressibility modulus when the monolayer is in the LC phase, showing that the monolayer is rigid at smaller A MA , before the monolayer collapses and the compressibility modulus drops again.</p><p>In experimental monolayer studies on MAs, the compressibility modulus indicates the liquid condensed (LC) phase with a value between 80-140 mN.m -1 , while a value between 20 -80 mN.m -1 indicates the liquid expanded (LE) phase. 29 The compressibility modulus for both systems simulated compares well with experimental values at large areas when the monolayers are in the LE phase. The compressibility modulus peaks at 133.15 mN.m -1 for the large MA monolayer while in the LC phase, falling within the range of experimental values. However, the small MA monolayer peaks outside of the experimental range at 238.59 mN.m -1 when in the LC phase. The exaggerated rigidity at smaller areas for this monolayer is likely caused by the larger percentage of MAs expelled from the water surface, forming an additional layer at the vacuum surface.</p><!><p>In experimental monolayer studies with MAs at large A MA , the molecules are proposed as folded at each functional group with four alkyl chains parallel to each other. As the monolayer is compressed, the molecule is proposed to extend its mero-chain, thereby decreasing its surface area to 3 or 2 parallel chains. A model for the unfolding of MA is shown in Figure 2.</p><p>It is possible to see the distribution of conformers relative to the structures in Figure 2 by creating a principal component (PC) map of the structures in Figure 2 A and then projecting the MAs from simulations onto the same space (Figure 9 C). On the PC plot, principal component 1 provides resolution on MA conformers in terms of the number of parallel chains. Lower to higher values on the axis move from the W-fold with four parallel chains to the Z-folds with three parallel chains, and finally the U-fold with two parallel chains prior to the fully extended single chain conformation. The sZ conformer has large overlap between the chain termini, and therefore has three parallel chains for most of its length, in contrast to eZ with its single chain section. Therefore, sZ is shifted more to the left relative to eZ. The relative energy of the conformers, shown in Figure 9B, indicates that the proximity of parallel chains leads to more stable conformers with lower energies. The W-fold with four parallel chains is lowest in energy, Z-conformers ranged from 10.52 -36.30 kJ.mol -1 , eU significantly higher at 89.48 kJ.mol -1 and the completely extended conformer has the highest relative energy of 94.13 kJ.mol -1 . On the second principal component axis the conformers are separated further, with an approximate zero value on this axis representing a conformer with no folding at any functional group, while negative values characterise conformers folded at the proximal cyclopropane (c) in contrast to conformers that do not fold at the group c, located at positive values on this axis. The PC plot shows that even at large surface areas the molecules did not stay in a four-chain folded conformation, which was the starting conformation. This finding contradicts current popular belief that MAs assume a folded conformation, such as a W-fold at large surface area. In fact, when taking the average structure of all the molecules in each surface area, the molecules are seen to be stretched out from the start, only bending at the acid group, which is in contact with the water surface. This finding that the AMA mero-chain is largely unfolded is in line with AMA being flexible under alkane conditions 24 and correlates with MAs folding more rigidly into W-conformations when an -methyl trans-cyclopropane group is present rather than cis-cylopropane. 32,33 At large areas of A MA = 1.02 -0.71 nm 2 the two chains are separated from each other, with most MAs having a conformation intermediate between eU and the extended straight conformer. The average structures here show a distance between a and d between 2.38 -3.64 nm. As the surface area is decreased, the chains come closer to each other. From A MA = 0.71 -0.61 nm 2 , the distribution of MA conformers on the PC plot starts to narrow down to the eU conformer, and at A MA = 0.51 -0.41 nm 2 there is a narrow distribution around the eU-conformer.</p><p>This shows that at A MA = 0.51 -0.41 nm 2 the MAs are tightly packed in the monolayer. As shown in Figure 10 for the large monolayer system, the order parameter for bonds in the alpha-chain region, C6D until C1D and the mero-chain region between the two cyclopropane beads M1B to M1A increases steadily from</p><p>A MA = 1.02 nm 2 where S ≈ 0 to A MA = 0.51 nm 2 where S ≈ 0.6. This indicates that the orientation of the bonds changes toward a parallel orientation with the monolayer normal as A MA decreases. From A MA = 0.41 nm 2 and smaller, the order parameter for these bonds decreases again until reaching a low of about 0.1, indicating random orientations at the end of the monolayer compressions. The point where S starts to decrease correlates with the collapse point. The order parameters for bonds in the short chain between the acid head group and the proximal cyclopropane (C3C-C1C) remain quite low with a maximum around S ≈ 0.2.</p><p>The bonds at the end of the mero-chain, from C5A to C1A increase very little in order between A MA = 1.02 -0.71 nm 2 , with S ≈ -0.2, indicating that the bonds are tending toward a perpendicular orientation, which means the end of the merochain is not stretched out vertically at these surface areas. Only from A MA = 0.61 nm 2 is a more significant increase in the order parameter seen. Again, a maximum is reached at A MA = 0.51 nm 2 before a rapid decrease. However, for this part of the chain the highest value reached is S ≈ 0.24, which is much lower than the order parameter for the alpha-chain and the proximal part of the merochain. Therefore, the parts of the alpha-and mero-chains that match each other in length seem to aid in increasing the order of those regions by supporting each other, while the longer terminus of the mero-chain is not matched and therefore does not pack very rigidly, even when the monolayer is highly compressed. This observation confirms previous interpretations from experimental monolayer studies that the length of the alpha-chain is important in how MAs fold, 27 as well as the relative length of the chain segments in relation to each other. 30,31 The order parameters for bonds at the functional groups in the molecule are lower than those for the C1-bead chains, except for the short chain between the acid head group and the proximal cyclopropane (C3C-C1C). S is very similar for M1B-C4B and C1B-M1A bonds, but C1C-M1B is generally higher than M1A-C5A.</p><p>Around the acid head group S stays between -0.2 and 0, indicating perpendicular or random orientations of the bonds in this segment.</p><p>The trends above are also seen in the small monolayer system. For the alpha and mero-chain segment, between M1A and M1B, S increases up to about 0.6 before decreasing again from A MA = 0.34 nm 2 , similar to the large system.</p><!><p>We present here a coarse-grained model for an AMA from M. tb.</p><p>Parameterisation of bonded interactions was performed on a single AMA in different environments, using atomistic MD as reference. The coarse-grained model replicated the atomistic model well. The peripheral OOH-COH-C1D</p><p>angle distribution was less precise, in favour of the angles that describe the alkyl chain conformation. The two cyclopropane beads were treated as identical, although their atomistic angle distributions varied slightly. This approximation resulted in the smaller angles for the distal cyclopropane bead, M1A being neglected. Hence, more folding may be expected at the distal cyclopropane group than is observed in this model.</p><p>Isotherms for the simulated homogeneous AMA monolayers were obtained for a small monolayer system containing 220 MAs as well as for a large monolayer containing 1972 MAs at 300 K. Despite the obvious difference with experimental monolayers that are heterogeneous in composition, even for single MA classes, the isotherms obtained from simulations were able to capture the key features of experimental monolayers such as phase transitions and monolayer collapse with reasonable accuracy. According to the calculated compressibility modulus, the monolayers were generally in the LE phase at higher A MA and in the LC phase at smaller A MA . A high compressibility modulus was obtained just before the monolayer collapse point, indicating a highly compressed monolayer, although a solid-condensed phase was not observed, in accordance with experimental data for AMA. 29,31 It was unexpected that MAs moved away from the water surface at larger A MA .</p><p>This observation may be due to an artefact in the model, or due to stabilising effects from the heterogeneity of experimental monolayers and the cell membrane. MAs moving away from the water surface resulted in systematic deviations in the A MA reported. These deviations were small at large A MA but increased as the monolayer was compressed, and more MAs moved away from the water interface. MAs were never seen to enter the water phase, unlike other lipids. 40,41 Investigation into MA conformations in the monolayer revealed that the postulated W-shaped folding at large A MA did not occur. Instead MAs were mostly folded only at the head group, with the two parallel chains becoming more tightly packed as the surface area decreased. It has been shown that the ciscyclopropane group does not promote folding to the extent that the -methyl trans-cyclopropane group does. 32,33 However, our model may be showing less folding at the cyclopropane beads than is occurring in reality, due to smaller angles being excluded in the CG model's angle distributions. Our model confirms the extension of the mero-chain as the surface area decreases and the molecules become more tightly packed, in accordance with experimental hypotheses. [25][26][27][28][29][30][31] The order parameter showed that tight packing with bonds orientated more perpendicular to the membrane normal was observed, as the monolayer was compressed, mostly for the alpha chain and the portion of the mero-chain adjacent and parallel to the alpha-chain. The mero-chain terminus, however, did not pack tightly, even at small A MA . This observation supports the theory that the alpha-chain length plays an important role in determining monolayer properties, as suggested by Hasegawa and co-workers. 27 To our knowledge this is the first coarse-grained model of a mycolic acid. The model cannot incorporate the explicit stereochemistry of a MA, which may be crucial in determining its structure and function. 8,9 Despite this, the model could replicate the key aspects of existing experimental data on MA monolayers.</p><p>Therefore, this approach can serve as a basis for further work on other MAs and their conformations, as well as models of larger components of the mycobacterial cell wall. All MA parameters are available for download on the Martini web portal cgmartini.nl.</p>