Template:Team:Groningen/CONTENT/modeling

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.<a href="#reps" class="cite">1</a> 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.<a href="#runw" class="cite">2</a> 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.<a href="#rdilsol" class="cite">3</a> 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 <a href="#pgachem" class="link">1</a>). 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.

<a href="#rbiomod" class="cite">4</a><a href="#rbiomod2" class="cite">5</a> 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.<a href="#rmartini" class="cite">6</a> 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.<a href="#rnafion" class="cite">7</a><a href="#rproton" class="cite">8</a>

An introduction to molecular dynamics can be found in the supporting material (<a href="Modeling/Supporting_material" class="link">link</a>), as well as some additional information on the creation of the protocols in the methods. In the methods section (<a href="#smethods" class="link">link</a>) 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).

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

0.28

12000

T3-M1-M2
M3-(+M1)-(+M2)
M3-T4-T5

114.646

25

M1

Nd

M1-M2

0.24

22000

T3-M1-M3
M3-(+M1)-(+M3)

124.643

250

M2

Qd

M1-M3

0.32

2000

M1-M3-(+M1)
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 <a href="#tparam" class="link">1</a>, while corresponding mapping can be seen in Figure <a href="#fmbeads" class="link">2</a>. 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 <a href="#fangles" class="link">3</a>. 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 <a href="#fe2e" class="link">4</a>, 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 <a href="#movaggr" class="link">1</a>). 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.<a href="#rwhc" class="cite">9</a> As can be seen in Figure <a href="#fcross" class="link">5</a> 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.

<video class="image" controls=""> <source src="https://static.igem.org/mediawiki/2015/3/32/Igem.groningen.2015.figure.small.modaggregation.mp4" type="video/mp4"> Your browser does not support the video tag. </video>

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.

<video class="image" controls=""> <source src="https://static.igem.org/mediawiki/2015/d/d9/Igem.groningen.2015.figure.small.modboxrot.mp4" type="video/mp4"> Your browser does not support the video tag. </video>

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 <a href="#movboxrot" class="link">2</a>. The concentrations in the simulations are equal to the concentrations used in the lab. In Movies <a href="#movpga" class="link">3</a>, <a href="#movcel" class="link">4</a> and <a href="#movcelp" class="link">5</a> 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.

<video class="image" controls=""> <source src="https://static.igem.org/mediawiki/2015/8/8d/Igem.groningen.2015.figure.small.modpga.mp4" type="video/mp4"> Your browser does not support the video tag. </video>

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.

<video class="image" controls=""> <source src="https://static.igem.org/mediawiki/2015/9/9f/Igem.groningen.2015.figure.small.modcel.mp4" type="video/mp4"> Your browser does not support the video tag. </video>

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.

<video class="image" controls=""> <source src="https://static.igem.org/mediawiki/2015/4/41/Igem.groningen.2015.figure.small.modcelp.mp4" type="video/mp4"> Your browser does not support the video tag. </video>

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 <a href="#fnaclw" class="link">6</a>. 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 (<a href="http://www.gromacs.org/" class="link">external website</a>, 2015) was used with verion 2.2 of the MARTINI coarse grained force field (<a href="http://md.chem.rug.nl/cgmartini/" class="link">external website</a>, 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.<a href="#rmartcarbo" class="cite">10</a> 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.<a href="#rmartcellfiber" class="cite">11</a> 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 <a href="#tcelparam" class="link">2</a>. The timestep used in all cellulose simulations was 20 fs.

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: <a href="Modeling/Supporting_material" class="link">link</a>

Marvasi,Massimiliano and Visscher, Pieter T. and Casillas Martinez, Lilliam

Exopolymeric substances (EPS) from Bacillus subtilis: polymers and genes encoding their synthesis

FEMS Microbiology Letters

1

1-9

2010

Stanley, Nicola R. and Lazazzera, Beth A.

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

Nagy, Z. and Novák, L. and Kozma, C. and Berka, M. and Bányai, I.

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

Horn, Harald and Lackner, Susanne

Modeling of Biofilm Systems: A Review

Productive Biofilms

ISBN: 978-3-319-09694-0

53-76

2014

Esser, Daniel S. and Leveau, Johan H.J. and Meyer, Katrin M.

Modeling microbial growth and dynamics

Applied Microbiology and Biotechnology

issn: 0175-7598

1-16

2015

Siewert J. Marrink,*,† and H. Jelger Risselada,† and Serge Yefimov,‡ and D. Peter Tieleman,§ and and Alex H. de Vries†

The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations

The Journal of Physical Chemistry B

27

7812-7824

2007

Aleksey Vishnyakov and and Alexander V. Neimark*

Molecular Dynamics Simulation of Microstructure and Molecular Mobilities in Swollen Nafion Membranes

The Journal of Physical Chemistry B

39

9586-9594

2001

E. Spohr,* and P. Commer, and and A. A. Kornyshev

Enhancing Proton Mobility in Polymer Electrolyte Membranes: Lessons from Molecular Dynamics Simulations

The Journal of Physical Chemistry B

41

10560-10569

2002

Na-Ri Lee and Tae-Hun Go and Sang-Mee Lee and Seong-Yun Jeong and Geun-Tae Park and Chang-Oh Hong and Hong-Joo Son

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

Cesar A. López and Andrzej J. Rzepiela and Alex H. de Vries and Lubbert Dijkhuizen and Philippe H. Hünenberger and Siewert J. Marrink

Martini Coarse-Grained Force Field: Extension to Carbohydrates

Journal of Chemical Theory and Computation

12

3195-3210

2009

César A. López and Giovanni Bellesia and Antonio Redondo and Paul Langan and Shishir P. S. Chundawat and Bruce E. Dale and Siewert J. Marrink and S. Gnanakaran

MARTINI Coarse-Grained Model for Crystalline Cellulose Microfibers

The Journal of Physical Chemistry B

2

465-473

2015