Team:Groningen/Modeling
Blue Bio Energy
Modeling poly-gamma-glutamic acid as a cation exchange membrane
Using classical molecular dynamics (MD), it is shown that poly-γ-glutamic acid can be used as a cation exchange membrane. This is done by comparing membranes made of different materials. Firstly, poly-γ-glutamic acid (PGA) was parametrized in the MARTINI Coarse Grained force field and our MD simulation shows that PGA molecules aggregate in dilute solutions. Consequently, we simulate reverse electrodialysis cells using aggregates of PGA, cellulose and cellulose phosphate, and we compare their performance. The simulation results show that the cellulose does not hinder the flow of ions and the cellulose phosphate acts as a strong cation exchange membrane, which conforms with literature. On the other hand, PGA performs as a cation exchange membrane.
The goal of the Bio Blue Energy project is to use bacteria to make a cation exchange membrane. For this to work the bacteria needs to have two properties. Firstly they need some way to stick together. If the bacteria can stick together, then a membrane can build from them. Secondly they need to filter for cations. Without this property the bacteria would just block or allow everything in the water to flow through and it would not be able to generate a current. Hence these are the crucial properties that need to be reproduced.
Some bacteria naturally have the first property. In this case they usually form a biofilm, which is an aggregate of bacteria in the slime they produce. One example of such a bacteria is the wild type Bacillus subtilis. While this bacteria makes a biofilm, our experiments show that this biofilm is not robust enough to withstand being in water for an extended period. In addition for the time that the biofilm is still intact it does not function as an ion exchange membrane.
Now, if one wants to add or improve the desirable properties to the biofilm, it is required to understand how the biofilm works. After all if one does not, then it is hard (if not impossible) to make a decent hypothesis to base the experiments in the lab on. What makes understanding the biofilm hard is that the exact composition of the the biofilm is not known, since it depends on the B. Subtilis strain, what genes are expressed, how the bacteria specialize and how heterogeneous the biofilm is.1 In addition, some genes related to the production of the biofilm are not well understood. While some are known to be responsible for the production of extracellular polymers, the length of the polymers in the biofilm varies a lot. This varying length makes characterization of these polymers hard.
While studying the literature on B. Subtilis biofilms we found that poly-γ-glutamic acid (PGA) might help us with both properties. It helps B. Subtilis stick together, because PGA helps with the underwater growth.2 This happens probably by increasing the number of sites available for salt bridges. It may also help the aggregation of B. Subtilis, because PGA forms aggregates even in dilute solutions.3 Poly-γ-glutamic acid may also form an obstacle for some ions, because is has a high number of fixed charges (as can been seen in its structure as shown in Figure 1). With a model of PGA one could gain insight in how it aggregates and how it could improve ion conductivity. Naturally there are multiple ways to approach modeling, usually though one takes a top-down approach or a bottom-up approach.
The structural formula of poly-γ-glutamic acid
When using the top-down approach, the models are usually based on the biomass, individuals (bacteria) or particles. These are useful to study the detachment, mass transfer and species distributions in biofilms, but fail when they are applied in bigger systems like whole reactors.
45 Furthermore, many parameters required for such models are hard to gather experimentally. Since this project is concerned with a complete system like a reverse electrodialysis cell and the electrostatic interactions, which are not included, this type of model would not be a good fit. While there are models for the growth of the biofilm, there are also models about the ion permeability of membranes. These models give good results with materials that have known geometries, but are less reliable for materials with pore sizes that are small relative to the Debye length or if the materials in question are heterogeneous.
Our model is based on the the bottom-up approach. Since poly-γ-glutamic acid was identified as a molecule that might have a positive contribution to our modified biofilm, the properties of this polymer were modeled. With molecular dynamics one is able to model the structure and dynamics of all kind of systems made of atoms and molecules. Thus molecular dynamics is able to simulate the interactions between the PGA molecules and to simulate the interaction of PGA with ions. Unfortunately, these simulations take a lot of time and computational power. Luckily it is possible to reduce the time needed by using the MARTINI coarse grained force field due to Marrink et al.6 While ion exchange membranes and proton exchange membranes have been studied using molecular dynamics before, this has never been done using MARTINI, and certainly not for poly-γ-glutamic acid.78
An introduction to molecular dynamics can be found in the supporting material (link), as well as some additional information on the creation of the protocols in the methods. In the methods section (link) protocols are detailed for the parameterize poly-γ-glutamic acid, the aggregation molecules, the membrane formation and the reverse electrodialysis cell simulations. In the results section we will report our findings and discuss them, and we will also mention possible future work. In the supporting material an overview of the specific tools and scripts we have written to create, modify or analyze topologies and trajectories is also given. And finally we give our conclusion.
Results
The modeling of poly-γ-glutamic acid consisted of the carrying out the following steps in order: parameterization of PGA, aggregation of PGA, membrane formation and ion flow simulations with PGA, cellulose and cellulose phosphate. Our main efforts in the parameterization of poly-γ-glutamic acid were focused on the bonded parameters. Interactions are bonded if they are part of a bond, an angle or a dihedral. Note that this property is not transitive.
The mapping of PGA to coarse grained beads. Each circle represents a particle (bead) in the coarse grained simulation. The atoms these circles encompasses in the structural formula of poly-γ-glutamic acid are part of the bead. The naming of these beads is based on if the bead is part of the monomer (M) or a terminal (T).
Bead
Type
Bond
Length
Force constant
Angle
Angle
Force constant
T1
Qd
T1-T2
0.23
18000
T2-T1-T3
111.013
100
T2
Qa
T1-T3
0.33
4000
T1-T3-M1
141.237
55
T3
Na
T3-M1
M3-(+M1)
M3-T4
M3-(+M1)
M3-T4
0.28
12000
T3-M1-M2
M3-(+M1)-(+M2)
M3-T4-T5
M3-(+M1)-(+M2)
M3-T4-T5
114.646
25
M1
Nd
M1-M2
0.24
22000
T3-M1-M3
M3-(+M1)-(+M3)
M3-(+M1)-(+M3)
124.643
250
M2
Qd
M1-M3
0.32
2000
M1-M3-(+M1)
M1-M3-T4
M1-M3-T4
133.794
250
M3
Na
T4-T5
0.24
4000
M2-M1-M3
113.393
40
T4
Nd
T4-T6
0.36
4000
M3-T4-T6
127.488
100
T5
Qa
T5-T4-T6
116.572
45
T6
Qa
The bonded interaction parameters of PGA. The Qd bead types have a single positive charge and the Qa types have a single negative charge. The beads with names starting with a "T" are part of the terminals of the polymer, while "M" denotes that it is part of the monomer. A plus sign in the name of bonds or angles denotes that it refers to a bead in the next monomer. The bond lengths are given in nm and the bond force constants in \(\text{kJ} \cdot \text{mol}^{-1} \cdot \text{nm}^{-2}\). The angles are given in degrees and the angle force constants are in \(\text{kJ} \cdot \text{mol}^{-1}\).
The final parameters for bonded interactions of PGA in the coarse grained model can be seen in Table 1, while corresponding mapping can be seen in Figure 2. The bead types were chosen based on expected hydrophobic/hydrophilic properties of the chemical groups. However changing the bead types of the backbone of the monomer did not have any effect on the end-to-end distance, bond and angle distributions(data not shown). In addition the shape of the bond distributions is similar but shifted, since the average bond lengths are a little bit higher in the coarse grained case (data not shown).
The angle distributions for beads in PGA terminals and monomers. The red graphs are from the all atom simulations and the green graphs are from the coarse grained simulations. Each subfigure shows an angle or a group of angles. These graphs show that the fit of the coarse grained model to all atom one is fine.
The distributions of end-to-end distances of PGA. When the parameterization is good, then the graphs have similar shapes. The end-to-end distance is a metric, which shows the distance between the two ends of a molecule. For polymers it is a measure of how open or closed the conformation is. Hence when the distributions of end-to-end distances are alike, then this means that the conformations of the molecules appear with similar frequencies.
The angle distributions are shown in Figure 3. The fit is good for the shape the all atom distributions have, it impossible to approximate the distribution perfectly in MARTINI. The reason for this is that the mapping of all atom data to beads may result in a non-harmonic potential, while the coarse grained angles are defined by a harmonic potential. Changing the force constants of the angles did not affect the distributions much, except if the constants were changed to a big value like 500 kJ/mol. However if the angle force constants are this big, then the distributions are too narrow. More importantly the distribution of the end-to-end distances, as can be seen in Figure 4, are comparable with each other. The only differences are that the the peak around 0.9 nm in the coarse grained graph is more distinct than in the all atom one. Furthermore, the distribution is narrower. The latter difference together with the angle distributions indicates that either the interactions between beads are too strong or the interactions of the backbone with the water is too weak. To fix this one would have to change the nonbonded interactions of the beads in PGA.
Membrane cross sections of A) Cellulose, B) PGA and C) Cellulose phosphate. Gray parts are the membrane and the blue parts are water. The cellulose membrane is denser than the other two, since the cellulose is in some parts of the membrane in a crystalline state. The cellulose phosphate is more more open than PGA, likely because of the bigger repulsion between the fixed charges.
Next, if multiple molecules are simulated together, they aggregate (as shown in Movie 1). As can be seen in the movie a shell of sodium ions forms around the polymers and some ions stay in position for a while. If the residency time is high enough, then this is counterion condensation. This phenomenon is usually observed if the coulomb interactions are stronger than the thermal interactions.
From an aggregate membranes were made and the water-holding capacity of the membrane made of poly-γ-glutamic acid was compared with values found in literature.9 As can be seen in Figure 5 the amount of water for PGA is approximately in agreement with the expected water-holding capacity of 56.9%. This property was not checked for cellulose phosphate membranes. Cellulose was already parameterized as a fiber, so validation was not needed. Comparing the other two cross sections with PGA one can see that the PGA membrane contains more water than one from cellulose, but less than one from cellulose phosphate. Note that the higher water content does not necessarily mean that the water and/or ion flow is less impeded.
The aggregation of PGA 5 molecules. Orange points are sodium ions, green points are chloride ions and the gray (neutral), purple (negative) and blue (positive) is part of PGA. 125 ns are shown in which all the molecules aggregate.
The simulation box of PGA. The membrane is gray, the sodium ions orange, the chloride ions green and the water blue. The salt concentration in the middle is 17 mMol (1 g/L) and the concentration in the outer sections is 510 mMol(30 g/L). The molecules and ions are shown in the initial condition of the ion flow simulation.
Finally simulations were prepared with two membranes consisting of a single type of molecule. A typical example of this kind of system can be seen in Movie 2. The concentrations in the simulations are equal to the concentrations used in the lab. In Movies 3, 4 and 5 the first 50 ns of the simulations are shown. These simulations show what kind of influence the fixed charges have on the ion conductivity. Based on these movies the sodium ions seem to be able to hop between the charge groups and move through the counterion shells in the membrane to the fresh water.
Ion flow through PGA membranes. The first 50 ns of the simulation are shown. Orange points are sodium ions, green points are chloride ions and the gray is PGA. As expected only sodium ions move to the fresh water.
Ion flow through cellulose membranes. The first 50 ns of the simulation are shown. Orange points are sodium ions, green points are chloride ions and the gray is cellulose. As is shown sodium and chloride ions move to the fresh water as was expected.
Ion flow through cellulose phosphate membranes. The first 50 ns of the simulation are shown. Orange points are sodium ions, green points are chloride ions and the gray is cellulose phosphate. As expected only sodium ions move to the fresh water.
While the effectiveness of PGA can be compared based on these videos, it is hard to make definitive statements about the flow rate of ions. Therefore we also calculated the net number of ions moving from salt water to fresh water and vice versa and then normalized the resulting value over the total number of the respective ions. The result can be seen in Figure 6. As is shown PGA and cellulose phosphate let only sodium ions through while blocking the chloride ions. The cellulose allows both sodium and chloride ions to pass through. Another noteworthy thing is that in the beginning the PGA lets through the sodium at the same rate as cellulose. It is unclear if this is so because of the different charge, pore size or width in comparison with cellulose phosphate. In the figure the water flow is also shown. The graphs show that the water moves from the fresh water to the salt water. The main reason for this effect is that the water has to fill the space left by the ions moving from salt to fresh water. Because cellulose does not filter much at all, it might look strange that the water flow is so low for cellulose. The reason is that the total number of ions is lower in this simulation and therefore the total of number of water molecules moved to the salt water is lower too.
The ion and water flow through the membranes over time. Green is PGA, red is cellulose and blue is cellulose phosphate. Note that virtually no chloride ions pass the membrane with PGA or cellulose phosphate. In addition the sodium flow is slower with cellulose phosphate than with PGA. The cellulose membrane does not filter at all, since sodium and chloride ions pass the membrane in equal ratio's.
Future work
Before these results can be taken as granted, some points have to be addressed. Firstly, since the electrostatic interactions are not well represented in MARTINI, the ion flow simulations should be validated by backmapping the system to an all atom model and simulating the result. This would give insight in how well the MARTINI CG force field works for this type of simulation, since the interactions between the ions and the membrane is more accurate in the all atom simulation.
Secondly, since the properties of the poly-γ-glutamic acid model do not fit the ones calculated from the all atom simulations and experimental data, it would be nice if the parameterization is improved. For example the non-bonded interaction parameters could be investigated, since they also play an important role in the molecular behavior and were neglected in our parameterization.
Thirdly, the influence of charges in the membrane is not clear, since different charge configurations were not investigated. It might be that only charges at the surface of the membrane are important, or maybe denser or stronger charges improve the flux and/or selectivity.
Fourthly, the results of aggregation were not thoroughly checked. The length of the simulated polymers is rather small in comparison with the length of polymers used in experiments, which could lead the simulated molecules to exhibit different behavior. For example materials can behave quite differently at a nanoscopic scale then at a mesoscopic scale, because at mesoscopic scale the properties of the material is mainly defined by the properties of the aggregate rather than on the properties of single molecules. Therefore the aggregate of PGA should be checked for the properties it exhibits at mesoscopic scale.
Lastly, while poly-γ-glutamic acid is a promising substance, others might also provide the same benefits. Based on the result that charges in the molecule are important, one could come up with a few other molecules to build a membrane with. Other interesting options would be polyacrylic acid or poly-β-aspartic acid, since these have a higher charge density than PGA. Or maybe one could try to make a AEM instead by using ε-polylysine, because its structure is similar to PGA but with positive functional groups.
Methods
Extended explanations for the chosen methodology can be found in the supplementary material. The molecular dynamics software GROMACS 4.6.7 (external website, 2015) was used with verion 2.2 of the MARTINI coarse grained force field (external website, 2015). Also note that with "all atom" topologies based on the united-atom force field, while "coarse grained" refers to topologies based on the MARTINI force field.
Molecule parameterization
The parameterization of poly-γ-glutamic acid is based on the general procedure to parameterize a new molecule in MARTINI. In essence the following steps are performed:
1 do an all atom simulation
2 choose a mapping to map the all atom molecule structure to the coarse grained one.
3 using this mapping, decide on and/or calculate the parameters for the coarse grained molecule.
4 run a coarse grained simulation with the new parameters.
5 compare some statistics to see if the coarse grained model fits the all atom simulation data.
6 repeat 2-5 until the fit is good enough.
The timestep used in these simulations was 40 fs in the case of the MARTINI CG force field and 2 fs for the all atom simulation. These simulations ran for 40 ns and 45 ns respectively. The correctness of the fit was based on the distributions of the angles, bonds and end-to-end distance.
Cellulose
The parameters for the monomer of cellulose is taken from the paper due López et al.10 Like they did for the hexamer of amylose, the topology of the monomer of cellulose was extended by adding bonds, angles and dihedrals. The angle over the backbone is based on the paper on cellulose fibers, while the other extensions are based on the similarities in the cellulose monomer.11 Note that the first paper models cellobiose and the second one crystalline cellulose. In our model the presence of both states are preferred. The parameter extensions are given in Table 2. The timestep used in all cellulose simulations was 20 fs.
Bond
Length
Force constant
Angle
Angle
Force constant
Dihedral
Angle
Force constant
B4-(+B2)
0.518
30000
B2-B4-(+B2)
168
50
(+B1)-(+B2)-B4-B5
30
5
B4-(+B2)-(+B4)
168
50
(+B1)-(+B2)-B4-B6
-150
5
(+B3)-(+B2)-B4-B5
-150
5
The extra cellulose bonded interaction parameters. The bond lengths are given in nm and the bond force constants in \(\text{kJ} \cdot \text{mol}^{-1} \cdot \text{nm}^{-2}\). The angles are given in degrees and the angle force constants are in \(\text{kJ} \cdot \text{mol}^{-1}\). The dihedral angles are given in degrees and the dihedral force constants are in \(\text{kJ} \cdot \text{mol}^{-1}\), while the multiplicity is 1 for all dihedrals.
Cellulose phosphate
The definition of cellulose phosphate is based on the definition of cellulose. The main difference is that the bead type of bead B1 and B5 are changed to Qa with a charge of -2. All simulations with cellulose phosphate use a timestep of 20 fs.
Aggregate formation
For aggregation a box size of 10 nm x 10 nm x 10nm is used, except for cellulose phosphate the box was 15 nm wide in the z axis direction. In such a box 35 molecules with 20 monomers (including terminals) are placed. The general procedure is as follows:
1 enlarge the box around 1 molecule or extract the aggregate from a previous simulation.
2 add more molecules and fill the rest of the box with water.
3 replace some of the water with ions and antifreeze based on the desired concentration.
4 run an energy minimization, relaxation, and a MD simulation to simulate the aggregation of the molecules.
5 repeat steps 1 till 4 until the target number of molecules is in the simulation box.
Membrane formation
The procedure to make the membranes is as follows:
1 take the end state of an aggregation simulation.
2 turn the water into polarizable water.
3 force the formation of the membrane in the xy plane by applying a force of 1, 5 and 10 kJ/mol for PGA and 1 and 5 kJ/mol for cellulose based membranes.
4 relax the membrane by decreasing the force constant in reverse order.
Ion flow tests
For this last phase a semiisotropic pressure coupling type is used instead an isotropic type.
The procedure is as follows:
1 double the box with the membrane in the z direction.
2 run the script to adjust the concentrations.
3 run an energy minimization while keeping the membrane in place by force.
4 run a relaxation simulation until the sodium just about reaches the surface of the membrane on the fresh water side. Again the membrane is kept in place by force.
5 run the production simulation while keeping only some points of the molecules in the membrane fixed.
A model for PGA in the MARTINI CG force field was presented and it was shown that PGA can be used as a cation exchange membrane. Our methodology and analysis was applied to different molecules (cellulose and cellulose phosphate) to gain further insight in the performance of PGA membranes. In addition some of the aggregation properties of PGA were reported. Together these results answer our original question, namely if it could positively contribute to the ion conductivity of biofilms.
The supporting material can be found here: link
Exopolymeric substances (EPS) from Bacillus subtilis: polymers and genes encoding their synthesis
FEMS Microbiology Letters
1
1-9
2010
10.1111/j.1574-6968.2010.02085.x
Defining the genetic differences between wild and domestic strains of Bacillus subtilis that affect poly-γ-dl-glutamic acid production and biofilm formation
Molecular Microbiology
4
1143-1158
2005
10.1111/j.1365-2958.2005.04746.x
NMR Study of Poly(γ-Glutamic Acid) and Partially Benzylated Poly(γ-Glutamic Acid): Nanoparticles in Solution
Colloids for Nano- and Biotechnology
ISBN: 978-3-540-85133-2
200-208
2008
10.1007/2882_2008_112
Modeling of Biofilm Systems: A Review
Productive Biofilms
ISBN: 978-3-319-09694-0
53-76
2014
10.1007/10_2014_275
Modeling microbial growth and dynamics
Applied Microbiology and Biotechnology
issn: 0175-7598
1-16
2015
10.1007/s00253-015-6877-6
The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations
The Journal of Physical Chemistry B
27
7812-7824
2007
10.1021/jp071097f
Molecular Dynamics Simulation of Microstructure and Molecular Mobilities in Swollen Nafion Membranes
The Journal of Physical Chemistry B
39
9586-9594
2001
10.1021/jp0102567
Enhancing Proton Mobility in Polymer Electrolyte Membranes: Lessons from Molecular Dynamics Simulations
The Journal of Physical Chemistry B
41
10560-10569
2002
10.1021/jp020209u
In vitro evaluation of new functional properties of poly-γ-glutamic acid produced by Bacillus subtilis D7
Saudi Journal of Biological Sciences
2
153-158
2014
10.1016/j.sjbs.2013.09.004
Martini Coarse-Grained Force Field: Extension to Carbohydrates
Journal of Chemical Theory and Computation
12
3195-3210
2009
10.1021/ct900313w
MARTINI Coarse-Grained Model for Crystalline Cellulose Microfibers
The Journal of Physical Chemistry B
2
465-473
2015
10.1021/jp5105938