Skip to main content

Exploring the organic–inorganic interface in biosilica: atomistic modeling of polyamine and silica precursors aggregation behavior


Diatoms are a significant group of algae displaying a sizeable morphological diversity, whose underlying structure arises from nanopatterned silica. Extensive experimental evidence suggests that a delicate interplay between various organic components and polysilicic acid plays a crucial role in biosilica mineralization. Thus, gaining insight into the properties of this organic–inorganic interface is of great interest in understanding the mechanisms controlling biosilica formation over different length scales. In this work, we use all-atom Molecular Dynamics simulations to investigate the aggregation behavior of polysilicic acid and silica nanoparticles in solution in the presence of protonated long-chain polyamines with a focus on the nature of the driving forces mediating the organic–inorganic aggregation process. Our results show that electrostatic forces between organic and inorganic species are the dominant interaction responsible for largely preserving the structural integrity of the organic–inorganic aggregates in solution. Thus, aggregates involving electrically neutral polysilicic acid are fully dissolved in an aqueous environment, since hydrogen bonding and van der Waals interactions turn out to be not strong enough to keep the aggregates together. Our main simulation results are in qualitative agreement with in vitro experiments, so that we expect they can contribute to shedding light on the initial stages of biosilica mineralization in diatoms.


Diatoms have enormous ecological relevance and display a diversity of patterns and structures ranging from the nano- to the millimeter scale. More than 250,000 species of diatoms have been identified, mainly based on their shell morphology. Apart from their interest in biology, diatom nanotechnology is emerging as a new interdisciplinary area, joining efforts of scientists working in biology, chemistry, biochemistry, physics, biotechnology, and materials science. Possible applications of diatoms in biophotonics, microfluidics, sensorics, controlled drug delivery and even in the development of unique computer architectures have been proposed, see Ref. [1] for a review. Intensive experimental investigations over the past years have identified a specialized intracellular compartment, the silica deposition vesicle (SDV) that is considered to be responsible for the control of biosilica precipitation and pattern formation (see Refs. [2, 3] for general reviews). Of great interest has been the discovery of a large variety of organic (both soluble and insoluble) components inside the SDV, whose role in biosilica formation has only been partly clarified. Although the identification of these silica-associated biomolecules seems not yet complete, the organic components identified so far include long-chain polyamines (LCPA), zwitterionic silaffins, silacidins, cingulins, and chitin [4,5,6,7,8,9,10,11,12,13,14,15,16,17]. In particular, LCPAs and silaffins have been shown to strongly influence in vitro silicic acid polycondensation as well as precipitation of silica nanoparticles with diameters covering a broad range of few tens to few hundreds of nanometers.

Clarifying the underlying mechanisms of biosilica morphogenesis in diatoms (and other species)—involving a diversity of length and time scales—will provide deeper insights into how biological information is transferred from the atomic scale to the macroscopic scale, and will thus allow suggesting specific strategies for possible nanotechnological applications. Still, a genuinely scale-bridging approach is missing, due to the enormous complexity of the problem. Instead, a variety of computational methodologies has been used to address various physico-chemical processes at different length scales, ranging from quantum mechanical methods [18,19,20], through (atomistic and coarse-grained) classical Molecular Dynamics simulations [21,22,23,24,25,26] up to models based on reaction–diffusion and phase-field approaches [27,28,29] and kinetic Monte Carlo simulations [30].

Concerning the potential role of some of the organic components, the currently existing general hypothesis is that they may, on the one hand, speed up or even catalyze silicic acid polycondensation [8] and, on the other hand, provide an organic template for biosilica precipitation. This second idea has been pioneered by the work of Sumper [31], see also Ref. [32]. Although the underlying physicochemical mechanisms are not fully understood, it is generally assumed that silica formation is accompanied by self-assembly processes involving a diversity of microscopic interactions: electrostatic, hydrophobic, van-der-Waals, and hydrogen bonding between the different biomolecules and between organic constituents and polysilicic acid. Without entering into further details, we refer the reader to review articles by Hildebrand and Lerch [33], Foo et al. [34], and Perry and Keeling-Tucker [17].

In the specific case of the long-chain polyamines, which will be the target of our computational studies, they were also found to mediate silica formation in vitro, see Ref. [2]. The LCPAs require, however, the presence of polyvalent anionic species to trigger silica deposition.

Experimental studies by Sumper et al. [35, 36] have revealed the dramatic influence phosphate anions have in controlling the size of silica nanoparticles formed in an LCPA containing solution. In fact, silica formation was only triggered if phosphates were present in the solution containing LCPA and silicic acid. An interesting observation in the same study was the presence of a lag phase in the silica precipitation experiments, which the authors associated with the time-lapse needed for mono silicic acid to polymerize spontaneously. This issue of the initial steps towards silicic acid polycondensation and the role played at this stage by the organic components is one central open question in biosilica formation in diatoms.

Subsequent investigations by Brunner et al. [37] using NMR and dynamic light scattering on polyallylamine (PAA) further demonstrated that silica precipitation took place only above a threshold in the phosphate ion concentration \(\sim \) 0.3 [Pi]/[r.u.], clearly indicating the crucial role of the phosphate ions. A similar correlation was found between the size of PAA aggregates and the ionic concentration with a dramatic increase in the diameter of the aggregates around the threshold value of 0.3. Interestingly, also sulfate anions were shown to induce phase-separation of PAA solutions, while chloride ions did not, even at very high ionic concentrations.

Lutz et al. [38] used NMR and dynamic light scattering methodologies to assess and characterize the phase separation behavior of PAA molecules in aqueous solution in the presence of phosphate ions with varying pH. This study revealed a strong pH-dependence of the average PAA aggregate diameter, with the threshold pH at which aggregation begins being dependent on the ion concentration. While these results indicated a strong influence of electrostatic interactions (controlled by the pH) in mediating the PAA self-organization, the authors also hinted at the presence of a hydrogen-bonded network stabilized by the electrostatic interactions. One must, however, keep in mind that PAA chains are considerably longer than LCPAs −LCPA molecules typically found in diatoms have a chain length ranging between 8 and 20 amino groups [8]—, so that the entropic contributions to the self-assembly may be quite different in both cases. Nevertheless, the relevance of electrostatic and hydrogen-bond interactions can be expected to be equally relevant for LCPA-phosphate ion interactions and hence, also control the morphology of biosilica.

The primary qualitative picture that seems to emerge from these results is that the positively charged polyamine chains interact electrostatically and eventually through hydrogen bonding with the polycations and with the (negatively charged) polysilicic acid molecules. This process results in the formation of silica networks composed of spherical aggregates with diameters depending on the specific reaction conditions.

The experimental studies previously reviewed build our motivation to carry out all-atom classical Molecular Dynamics simulations of biologically relevant organic components (LCPAs) in solution in the presence of phosphate ions and silicic acid nanoclusters (with diameters \(\le \) 1 nm) as well as amorphous silica nanoparticles (with diameters \(\sim \) 1–2 nm). Our main goal is to shed light on the possible self-aggregation pathways of organic and inorganic components, the aggregate’s structural stability in solution, and the dominant non-covalent interactions (electrostatic, hydrogen bond, van-der-Waals) controlling it. The choice of LCPAs is additionally motivated by detailed solid-state NMR spectroscopy investigations by Jantschke et al. [15] on S. turris, which strongly suggests that LCPAs are broadly dispersed into the silica phase. Moreover, the mass spectral characterization of polyamines from an Antarctic diatomaceous ooze sediment sample carried out by Bridoux and Ingalls [39] also indicated that LCPAs are tightly associated with diatoms cell walls. Hence, LCPA provides a relevant molecular system to explore the properties of the organic–inorganic interface. LCPAs are also very suitable for a computational study due to their relative structural simplicity, in contrast, e.g., to the zwitterionic silaffins, which possess a variety of post-translational modifications, making the modeling of large molecular assemblies computationally more expensive. As the bottom line of our computational study, our results strongly hint at the dominant role played by electrostatic interactions in influencing the structural integrity of LCPA/silica aggregates. Aggregates involving the electrically neutral dimers and tetramers are fully dissolved in an aqueous environment. We also discuss the influence of phosphate ions in influencing the self-organization of LCPAs in solution.


Object of study

Due to the intrinsic restrictions on the accessible time scales of all-atom Molecular Dynamics (MD) simulations, we aim in this study at revealing clustering tendencies—rather than a large scale self-assembly—of specific organic components in solution in the presence of inorganic components. We further aim at exploring the relevant physical interactions (electrostatic, hydrophobic, hydrogen bond, van-der-Waals) mediating the organic–inorganic clustering.

Organic components

As previously mentioned, we have chosen long-chain polyamines (LCPA) as the organic component due to both their relatively simple structure as well as to their relevance for biosilica formation. Figure 1a shows the atomistic representation of the LCPA to be used in this study. It contains eleven amino groups (a.g.), and it has been previously synthesized and characterized by Nuclear Magnetic Resonance experiments by Abacilar et al. [40].

The obtained pKa at physiological pH (\(\sim \) 5.5) indicates that the LCPAs are positively charged with a degree of protonation of \(9\, e\).

Fig. 1

Main components of the simulation a Long-chain polyamine (LCPA) containing 11 amino groups with a protonation degree of \(+9\). The two uncharged amino groups are highlighted with a blue triangle; b neutral disilicic acid (dimer|\({{\hbox {Si}}_2{\hbox {O}}_7{\hbox {H}}_6}\)), neutral cyclo tetrasilicic acid (tetramer|\({{\hbox {Si}}_4{\hbox {O}}_{12}{\hbox {H}}_8}\)), and negatively charged octahydroxyoctasilsequioxane-ion (octamer|\({{\hbox {Si}}_8{\hbox {O}}_{20}{\hbox {H}}_7^-}\)); c randomly generated amorphous silica nanoparticles


As mentioned in the introduction, in vitro experimental studies [41,42,43,44] have shown that phosphate anions have a dramatic influence on the size of silica nanoparticles that are formed in an LCPA containing solution. Since silica formation is only triggered if phosphates (or other polyvalent anions) are present in the solution, we have included phosphoric acid in our simulations. Since it is known [45] that more than 95 % of the species are present in solution in the form of dihydrogen phosphate ion (\({{\hbox {H}}_2\hbox {PO}_4^-}\)) at a pH of 5.5, we have used this charge state throughout this paper.

Silica precursors

We have selected different silica precursors with the purpose of shedding some light on the role of microscopic interactions at the organic–inorganic interface. The chosen silica precursors are displayed in Fig. 1b, c. In the first part of this study, we have used three different polysilicic acids as precursors: disilicic acid (\({{\hbox {Si}}_2{\hbox {O}}_7{\hbox {H}}_6}\)), cyclo tetrasilicic acid (\({{\hbox {Si}}_4{\hbox {O}}_{12}{\hbox {H}}_8}\)) and octahydroxyoctasilsequioxane-ion (\({{\hbox {Si}}_8{\hbox {O}}_{20}{\hbox {H}}_7^-}\)); for simplification, they are denoted in the text as dimer, tetramer, and octamer, respectively. Concerning the tetramer and the octamer, the cyclic structures have been considered because they turn out to be formed earlier in the condensation process [46] compared to the linear conformations. As we have verified with Marvin sketch, dimers and tetramers are neutral while the octamers are negatively (\(-1\,e\)) charged at pH \(\sim \) 5.5.

In the second part of the investigation, we have used larger amorphous silica nanoparticles to address the interactions with LCPAs on a larger length scale. As initial bulk material, the dataset generated by Cruz-Chu et al. [47] and published along VMD in the “Inorganic Builder” plugin was used. This set of coordinates was generated following the cooling cycles proposed by Huff et al. [48]. From the bulk material, random spheres with preselected diameters were extracted. Due to the nature of amorphous \({{\hbox {SiO}}_2}\), a triple coordinated oxygen atom may exist inside a randomly selected sphere. If this is the case, the sphere was discarded, as the used force fields can not describe the behavior of such a configuration. After a valid sphere is cut from the bulk, the system was “polished.” During this process, all silicon atoms having less than two connections to other silicon atoms were removed to create as compact a structure as possible. Next, all outward-facing oxygen atoms with just one connection to a silicon atom were identified. These selected atoms were then saturated with a hydrogen atom. In a first approximation, they were placed randomly with a distance of 0.95 Å and at an angle (\({{\hbox {Si}}{-}{\hbox {O}}{-}{\hbox {H}}}\)) of \({130}^{\circ }\). A random selection of hydrogen atoms was removed, and the connected oxygen atoms were charged to adjust the virtual pH value. Based on the pH dependence of the silica surface described in the INTERFACE force field 1.5 [49], we chose to remove 7.75 % of the hydrogen atoms to be near our pH \(\sim \) 5.5 goal. For the particles used, the charges ranged depending on the surface area from \(-1\) e to \(-4\) e. Charges in the structure were then set depending on their neighbors following the INTERFACE force field. After the “polishing” step the spheres were relaxed in vacuum at 298 K for 20,000 steps (2 ps).

Simulation and visualization details

The classical all-atom Molecular Dynamics (MD) simulations of LCPAs in solution with varying silica precursors and phosphate ion concentrations were performed using the software package NAMD 2.10.8 [50].

Polyamines and chlorine counter-ions were modeled using the CHARMM force field [51], while the silica species and phosphate ions were described using the INTERFACE force field 1.5 parameters [52]. The LCPA structure was first optimized in vacuum using Density-Functional Theory as implemented in Gamess [53] at the B3LYP/6-31G(d,p) level of theory. The atomic partial charges of LCPA were derived using the R.E.D (RESP ESP charge derive) server [54, 55]. Water molecules were described using the TIP3P model (with \(-0.834\) 4 charge for oxygen and 0.417 e for hydrogen). Bonds involving hydrogen were kept rigid during the simulations using the SHAKE algorithm [56]. Non-bonding interactions were computed for atoms separated by at least three bonds using a cutoff at 12 Å with the SWITCH [50] function turned off and using a scaling factor of 0.8 to reduce the electrostatic forces. Finally, long-range electrostatic interactions were evaluated using the particle-mesh Ewald [57] method, as implemented in NAMD with a grid spacing of 1.0 Å.

All the MD simulations were first equilibrated in an NPT ensemble for 5 ns to 15 ns. Once a reasonable value of the density is reached (0.9 g cm\(^{-3}\) ), consistent for organic components in the diluted water solution, the volume was kept constant carrying on the simulations in the canonical ensemble (NVT). All simulations were performed with an integration time step of 2 fs at temperature \(T={298}\,\mathrm{K}\) and pressure \(P={1.013}\, {\mathrm{bar}}\). The temperature was kept constant using a Langevin thermostat [58]. The simulation boxes are cubic and rectangular cuboids with periodic boundary conditions applied in all directions. The trajectories were collected for at least 50 ns saving every 1000 to 3000 steps. For the visualization of the molecular structures and their analysis, we used VMD [59] to generate the scenes that were then rendered with Tachyon [60]. To generate the animations (Additional files 1, 2, 3 and 4) the trajectories were prepared with the help of MDAnalysis [61, 62]. The system was reduced to the relevant components, and the building blocks paths were smoothed to capture the dominant conformational changes better. In the end, the rendering was carried out automatically with Molani [63].

Results and discussion

Interaction of polyamines with polysilicic acid

We have first addressed the issue of the stability of LCPA clustering in the presence of polysilicic acid (dimers, tetramers, octamers) and phosphate ions, identifying the role played by each element and the main driving forces in the mutual rearrangement among the organic–inorganic components.

As Brunner and coworkers have shown [64], multivalent anions such as phosphates can efficiently induce microscopic phase separation in the silica precipitation in the presence of polyallylamine. As hypothesized, phosphate ions should facilitate the formation of larger molecular aggregates due to the expected attractive electrostatic interaction between positively charged LCPAs and anionic species together with hydrogen bonds network formation. We have, therefore, considered two different simulation setups: (1) LCPA aggregation in solution with varying phosphate ion concentration, and (2) LCPA aggregation in the presence of polysilicic acid species at a fixed phosphate ion concentration.

1. Aggregation of polyamines in solution with varying phosphate ion concentration

The simulation cell contained 20 charged (protonation degree +9) LCPAs with different phosphate ion (Pi\(^-\)) concentrations. We have considered three Pi\(^-\) concentrations per charged amino group, [Pi]/[a.g]\(^+\): 0.5, 1.0 and 2.0. Representative snapshots are reported in Fig. 2. As a reference, a system without phosphate ions, including only LCPAs in solution, was also studied.

To reduce the computational effort, we have adopted the following strategy: before addressing the system in an aqueous environment, we have first pre-aggregated the organic and inorganic components in a vacuum; only then the obtained cluster was introduced in the simulation box including the water molecules and its structural stability tested as a function of the simulation time. In this way, we overcome the very long time window that would be required to induce clustering in solution starting from randomly placed individual components (in the absence of water, aggregation takes place very quickly due to the non-screening of Coulomb interactions). To avoid spurious results related to specific clustering configurations in vacuum, we repeated this procedure for different initial configurations, but the obtained final structures displayed similar structural features.

Fig. 2

Representative snapshots of the aggregation of 20 protonated LCPAs (green) at different phosphate anion (Pi\(^-\)) concentrations (red). The phosphate ion concentrations are: [Pi]/[a.g.]\(^+= 0.5\) (a, d); 1.0 (b, e), and 2.0 (c, f). The label “initial” (panels ac) refers to the starting pre-aggregated system in vacuum for each corresponding Pi\(^-\) concentration. The aggregate is subsequently immersed in an aqueous solution to test its structural stability. For clarity, the water molecules are not shown. Notice that only for the cases e and f the LCPA-phosphate aggregate remains partially intact. On the contrary, in case d the aggregate completely dissolves after a short simulation time (well below the shown snapshot at t = 20 ns)

In Fig. 2 we show initial (a–c) and final (d–f) configurations of the polyamine-phosphate clusters for increasing ion concentrations: [Pi]/[a.g.]\(^+= 0.5\) (a, d); 1.0 (b, e), and 2.0 (c, f). For panel (d) the final snapshot was taken after \(t=20{\hbox { ns}}\), while for panels (e) and (f) the final conformation was drawn at time \(t=70{\hbox { ns}}\). From the figures, we see that for the smallest concentration [Pi]/[a.g.]\(^+= 0.5,\) the pre-aggregated cluster is structurally unstable upon immersion in the water box and rapidly dissolves. However, the process slows down with increasing phosphate concentration. In fact, after 70 ns (see Fig. 2e, f), few LCPAs are dispersed in water, but nevertheless an internal core still remains intact for the entire simulation time. By analyzing in more detail the obtained simulation trajectories, we found out that the positively charged amino groups in the polyamines ionically cross-link with the negatively charged oxygens belonging to Pi\(^-\). Furthermore, a hydrogen bond network among phosphate groups additionally contributes to stabilizing the cluster, see Fig. 3. These results are in agreement with experimental studies [44, 65], which suggest not only Coulomb interactions as a stabilizing agent of the polyamine-phosphate network, but also hydrogen-bond interactions.

Notice that for [Pi]/[a.g.]\(^+\)= 2.0 the electrostatic repulsion effect among the phosphate ions may be too strong (a larger number of ions are dispersed in the solution) and, hence, the structural stability of the aggregate might become compromised over longer time scales not reached by the simulations. In this sense, the optimal concentration seems to be [Pi]/[a.g.]\(^+\)= 1.0. This effect may, however, depend on the concentration of polyamines involved in the simulation.

Fig. 3

Hydrogen bond network between LCPA and phosphate ions. The interaction mechanism between LCPA and Phosphate ions is shown in detail; specifically, the positively charged amino groups ionically cross-link with the negatively charged oxygens of Pi\(^-\) (see black dotted lines). In addition, the hydrogen bond network among phosphates is also highlighted (red dotted lines)

The observed trend with increasing phosphate concentration also agrees with the experimental studies in Refs. [8, 42, 44, 64, 66, 67], which demonstrates a delicate interplay between phosphate concentration and the solution pH (which allows tuning the charge state) in controlling the aggregation behavior of organic components and the subsequent silica deposition process. We also remark that, besides the biosilica-relevant experimental studies mentioned above, phosphate-induced self-assembly of polyamines into supramolecular nano-architectures has also been addressed to develop strategies for novel functional materials [65, 68]. These studies have also highlighted the role of electrostatic and hydrogen-bonding interactions with the phosphate ions in the formation of the polyamine molecular networks.

2. Aggregation behavior of polyamines in the presence of polysilicic acid at fixed phosphate ion concentration

In the previous section, we have shown the increased structural stability of the aggregates with increasing phosphate ion concentration. Starting now from the concentration [Pi]/[a.g.]\(^+=1\) for which we first observed a tendency to at least partly preserve the structural integrity of the polyamine/phosphate clusters, we have now considered three different systems corresponding to the three different polysilicic acid molecules shown in Fig. 1b.

As an initial configuration before immersion in the water box we have used in each case two different variants: (i) adding polysilicic acid to a pre-aggregated polyamine/phosphate system in a vacuum (as described in the previous subsection) and (ii) adding polysilicic acid to a randomly distributed set of polyamine and phosphate ions in a vacuum. In this way, we aim at clarifying if, and how, the initial conformation may bias the behavior in water solution. In both cases (i) and (ii) an MD simulation was subsequently run in vacuum to obtain an aggregate consisting now of polysilicic acid/polyamine/phosphates. This structure served then as the initial configuration for the simulations in an aqueous environment. These different stages are shown in Fig. 4.

Fig. 4

Two ways to generate the initial configuration of the polyamine (green)-silica dimers (grey)-phosphate (red) aggregate previous to the immersion in the water box. a Silica dimers are added in vacuum to a previously aggregated LCPA-phosphate complex. After an MD run, a compact cluster containing the three components is quickly formed, b cluster after immersion in water at 0 ns (initial configuration), c Snapshot showing the final stage of the MD simulation in water at 40 ns. The silica dimers and several phosphate ions are largely dispersed in the aqueous solution, while a polyamine-phosphate core is not fully dissolved. d All building blocks are randomly distributed in the simulation box without water molecules, e similar to the pre-aggregated case, a compact cluster is quickly formed and used as starting configuration for the MD simulations in water, f after 40 ns the aggregate is almost entirely dissolved in the solution

The transition from the structures shown in panel Fig. 4a, d to those in panel Fig. 4b, e (simulations in vacuum) took place just within a few nanoseconds (\(< {10}{\hbox { ns}}\)). Once the new aggregates are formed, we study their stability in water at ambient conditions.

Let us start by considering the case of the silica dimers (see Fig. 4) in the case of pre-aggregated initial conditions (panels a–c). Upon relaxation in a vacuum, the dimers distribute around a central core consisting of phosphate ions and polyamines (Fig. 4b). Leaving the system free to evolve in water leads to the dimers quickly dispersing (\(< {20} {{\hbox { ns}}}\)), see Fig. 4c, while the internal pre-stabilized polyamine/phosphate core remains relatively stable during the simulation time (60 ns).

In the case of a random initial distribution (panels d–f), all the components are obviously stronger mixed, with the dimers now being distributed inside and outside the final cluster (see Fig. 4e). Still, this initial conformation dissolves entirely once introduced in the water box (see Fig. 4f). This dissolution process is probably due to the presence of the neutral dimers, which screen the electrostatic interaction between polyamine and phosphate ions. These results suggest that the aggregation process is primarily driven by the Coulomb interaction between positively charged amino groups and negatively charged phosphate ions without any sizeable contribution of hydrogen bond among the neutral species. The same tendency was found in the presence of neutral tetrasilicic acid (not shown). This tendency may be an indication that the classical force fields used are not able to catch some subtleties of the interactions between neutral species and organic components. If this is the case, either more elaborate force field parametrizations will be required, or alternatively, a QM/MM approach may turn out to be more appropriate. This will require, in any case, further elucidation.

Let us now consider the system in the presence of negatively charged octamers at physiological pH \(\sim \) 5.5. Also in this case, the two previously described different pathways to generate the initial aggregates have been considered: pre-aggregated and random distribution, see Fig. 5.

Fig. 5

Two ways to generate the initial configuration of the polyamine (green)-silica octamer (grey)-phosphate (red) aggregate previous to the immersion in the water box. a Silica octamers are added to a previously built polyamine-phosphate cluster in a vacuum, b the resulting aggregate builds the initial structure at t = 0 ns for the MD simulation in a water box, c Snapshot of the final configuration after 40 ns. Although the phosphate ions are dispersed in the aqueous environment, an aggregate of polyamines and octamers remains structurally stable. df Similar results for the case that an initial random distribution of the three building blocks (polyamines, silica octamers, and phosphate ions) is used in a vacuum (d) to create the initial configuration for the simulations in water (e). After 40 ns, the phosphate ions are also dispersed in water as for c, but also, in this case, a mostly intact polyamine-octamer aggregate remains, in clear contrast to the simulations with silica dimers, see Fig. 4

The clusters containing octamers are formed in vacuum within few nanoseconds MD simulation, see Fig. 5b, e. In the case of the pre-aggregated initial condition, the silica precursors are distributed around the core of polyamine and phosphate ions. While in the case of random distribution, the aggregate is a mix of all components. Once the clusters are formed, we can again test their stability in solution. Representative snapshots are reported in Fig. 5. Independently of how the organic–inorganic aggregates were produced, we now see that even if some octamers are solvated, a significant portion of them are still part of the aggregate, which remains stable during the simulation time (40 ns), see Fig. 5c, f. Interestingly, for both initial configurations, a large number of phosphate ions quickly leave the aggregate, while the LCPA-octamer core shows a clear tendency to remain stable.

In order to characterize the obtained aggregates and determine the dominant interactions, the average number of contacts between LCPAs and polysilicic acid/phosphate ions have been computed using the Gromacs utility g-mindist [70] that calculates the number of contacts among molecules within a specific cut-off (0.6 ns). In our case, the number of contacts is divided by the total number of silica (180) or phosphate ions (180) and LCPAs (20). The results are displayed in panels a) and b) of Fig. 6 for silica dimers and octamers, respectively. As shown in Fig. 6a, the number of contacts between LCPA and silica dimers decreases quickly to zero in both configurations (see blue solid and dashed lines). This result confirms that interactions between LCPA and neutral silica components are not strong enough to keep the cluster together during the simulation time (see also the movie Additional file 1). In contrast, the number of contacts between LCPA and phosphate ions (green solid and dashed lines) is the only one significantly different from zero, resulting in the pre-aggregated polyamine/phosphate cluster remaining relatively stable during the simulation and the silica dimers quickly dispersing in the solution as we have previously mentioned.

Fig. 6

Time evolution of the average number of contacts between LCPA and either phosphate ions or silica precursors. The solid lines always refer to the random initial configuration while the dashed lines are related to the pre-aggregate initial configuration, see also Figs. 4 and 5 for reference. Contacts with phosphate ions are shown as green lines. a Dimers and b Octamers are used as silica precursors in separate simulations. The data was smoothed with a Savitzky–Golay filter [69]. The rapid decay in the average number of contacts between silica dimers and polyamines (a) reflects the rapid dissolution of the aggregate in the aqueous environment. The weaker decay in the average number of LCPA-phosphate contacts in a has most likely its origin in their mutual electrostatic interactions, which partly slows down the aggregate’s dissolution. In contrast, for the case of octamers (b), the average number of contacts between polyamine and silica remains nearly constant over the whole simulation time, illustrating the structural stability of the aggregates. However, the average number of contacts with the phosphate ions also decays as a function of the simulation time, but with a considerably slower rate

Concerning now the case of silica octamers, the situation is considerably different. In Fig. 6b we see that the number of contacts between LCPA and silica remains significantly different from zero during the entire simulation time (violet solid and dashed lines), while the number of LCPA/phosphate contacts tends to be reduced (green solid and dashed lines), although the decay is weaker than for the dimer case. More specifically, for pre-aggregate initial conditions (see the dashed lines in Fig. 6b), the number of contact increases in the case of LCPA/octamer (dashed violet line), but decreases for LCPA/phosphate (dashed green line). This trend suggests that some octamers replace phosphate ions. In the case of a random initial distribution, the number of LCPA/octamer contacts remains stable (solid violet line) along the simulation, while phosphate ions leave the cluster (solid green line). This suggests that once LCPAs and octamers are in contact, the configuration remains stable, and we do not find any exchange between octamers and phosphate ions.

These results show that van-der-Waals and hydrogen bond interactions between the silica precursors and LCPAs are in themselves not enough to stabilize the aggregates. Neutral precursors escape from the clusters and get solvated. On the other hand, in the presence of the charged silica octamers, the aggregates show a tendency to be stabilized, at least during the simulation time, suggesting that Coulomb interactions play a primary role in mediating the cluster stability. Our results are in agreement with a recent work by Centi et al. [26] claiming that the “neutral templating route” does not describe the early stage of the biosilification mechanism, but the ionic interactions among amines and silicates need to be considered.

As mentioned in the Introduction, an interesting observation in the experimental studies by Sumper et al. [42] was the presence of a lag phase in the silica precipitation process, which the authors associated with the time-lapse needed for mono silicic acid to polymerize spontaneously. A possible reason may be related to the fact that polysilicic acid becomes negatively charged only beyond a specific size (in our case from octamers on), and this will enhance the possibility for strong electrostatic interactions with the protonated LCPA molecules. Our simulation results presented so far point in a similar direction: only aggregates containing octamers show a tendency to remain structurally stable in solution. This is, however, not conclusive, since additional chemical interactions not accessible through classical force fields may play a role at the level of the charge-neutral species (dimers to tetramers). The role played by these interactions requires, therefore, a more in-depth inquiry. In particular, the hypothesized catalytic effect of polyamines [8] in the early stages of silicic acid polycondensation, facilitating siloxane-bond formation between monosilicic acid molecules, may play an important role and can only be addressed within a quantum mechanical framework.

Clustering of amorphous silica nanoparticles mediated by polyamines

From the previous discussion, it has turned out that, within our computational framework, the electrostatic interactions are by far playing the dominant role in determining the degree of structural integrity of the polyamine/phosphate ions/polysilicic acid aggregates. This result neglects the potential role chemical reactions play for the overall growth process. Newly formed chemical bonds could stabilize even small clusters formed just for a short time. As we can not address these interactions with standard classical force fields, we will focus in this section on clustering on longer length scales, involving larger silica nanoparticles. Here, we will assume that smaller precursors arising from aggregation processes involving polysilicic acid (on the length scales discussed in the previous section) have already formed tiny amorphous \({\hbox {SiO}}_2\) nanoparticles. The relatively large surface area of these particles results in multiple negative charges at the relevant pH \(\sim \) 5.5. Based on these assumptions, we investigate in the next step, the interaction of already preformed amorphous \({\hbox {SiO}}_2\) nanoparticles with LCPAs.

Since both components in this setup are multiply charged, the influence of the phosphate ions is less significant than during the early stages of aggregation. For this reason, we will not address the influence of phosphate ions, but only add a few additional chlorine ions to the system to achieve overall charge neutrality in the simulation box.

From the results presented so far in this study, we know that the preformed clusters consisting of negatively charged octamers, phosphate ions, and LCPAs can be relatively stable in water. To ensure that clustering directly in water is also possible, we started with a basic setup. In a water filled box (\(7 \times 6 \times 4\) nm) two spherical \({\hbox {SiO}}_2\) particles (NP1/NP2, \(d\approx {17}\) Å) are placed with a distance of 3.5 nm between their centers of mass. Between the two particles, a single LCPA is placed, as shown in Fig. 7a. At the desired pH 5.5, the LCPA has nine positive charged nitrogen atoms, and at the surfaces of the \({{\hbox {SiO}}_2}\) nanoparticles are six negative charges. Three chloride ions were added to make the box overall neutral. The overall system consists of 16,992 atoms that we simulated for 20 ns.

Fig. 7

Building blocks used for the simulation of the interaction between a polyamine chain and two silica nanoparticles. a Two \({\hbox {SiO}}_2\) nanoparticles (\(d\approx {13}\) Å) are placed around one LCPA molecule in a water box (\(7 \times 6 \times 4\) nm). b Snapshot of the emerging aggregate configuration after a simulation time of 20 ns

In the simulation (see Additional file 3 and Fig. 7), the LCPA is quickly attracted to one of the particles. After 2.5 ns, the molecule is wholly wrapped around one of the nanoparticles. The interaction is mediated by a few electrostatic anchor points and van-der-Waals interactions. The further dynamical evolution can be quantified by two quantities, which are shown in Fig. 8, where we plot the average distance between the centers of mass of the silica nanoparticles and the average distance of all atoms in the LCPA to the corresponding nanoparticle center of mass. Introducing these two quantities is necessary due to the structural flexibility of the LCPA. For the next 6.5 ns, no changes are seen until NP1 is, for the first time, sufficiently close to the NP2/LCPA group. At this point, the electrostatic forces are strong enough to affect the interaction between all three components. This interaction is possible since the charges on the LCPA are only partially compensated by the first \({{\hbox {SiO}}_2}\) particle. Between 14 and 16 ns, the molecule shifts more towards the NP1, and the system remains in this configuration for the remaining part of the simulation. In this basic setup, we observe fast clustering as soon as the components of the system are brought close enough together. Therefore, diffusion was the limiting factor in this simulation.

Fig. 8

Time evolution of the mutual separation of the three components. The distances between the two \({{\hbox {SiO}}_2}\) particles (NP1, NP2) as well as between each nanoparticle and the LCPA are shown during a 20 ns simulation run. The center of mass is used to measure the distance between the nanoparticles. For the LCPA, we use the average distance of its atoms to the center of mass of the corresponding nanoparticle. The near convergence of the three distances after 15 ns reflects the formation of the aggregate shown in Fig. 7

After we have established that clustering is possible, the interaction of multiple LCPAs with silica nanoparticles in the same system was investigated. To increase the number of components in our simulations while using a reasonable amount of computational resources was just possible by using a very high concentration of the raw elements. This high concentration reduces the distance the building blocks need to overcome before meaningful interactions can occur, therefore reducing the necessary simulation time significantly. The smaller box size also resulted in a lower number of water molecules.

In the first setup of this kind, 16 LCPAs are mixed with a variety of 24 \({{\hbox {SiO}}_2}\) particles (\(d\approx {10}\)– 20 Å). These building blocks are immersed in water in a cube with an edge length of 8 nm (see Fig. 9a). The charge difference of 89 e was compensated with chloride ions. Overall the system contains 49 606 atoms and the simulation was run for 62.8 ns. In a simulation with a large number of silica nanoparticles, an average distance loses its significance to represent the system. Instead, in the following, we track the size evolution of the newly formed clusters. A cluster is defined by a cutoff of 2.8 Å between the center of two adjoined atoms of different building blocks.

Fig. 9

Aggregation behavior of a polyamine-silica nanoparticle ensemble. a Snapshot of the initial configuration of 24 \({{\hbox {SiO}}_2}\) nanoparticles (gray; \(d\approx 10{-}{20}\) Å) with 16 LCPA molecules (green) in a 8 nm cube of water. b Snapshot of the final stage of the system after 63 ns. The supercell used is highlighted with a white square surrounded by repeating neighboring supercells to visualize the periodic boundary conditions used in the MD simulations

For every analyzed frame, the clusters were identified by extending next-neighbor searches. In Fig. 10 the results of this simulation are displayed. After the first 5 ns, nearly all available atoms are incorporated in a newly formed cluster. This fact remains unchanged for the rest of the simulation. The time evolution of the largest cluster in the system shows more detail in comparison. In the beginning, many small clusters exist that grow steadily by adding smaller clusters with just a few building blocks. As time goes on, nearly all available atoms are integrated into two to three large clusters. With a more substantial size, the successful combination of two clusters becomes more challenging. This temporary convergence of two large clusters results in sharp changes in the size of the largest one. An additional hurdle is the reduced mobility of the cluster due to the larger mass.

Fig. 10

Ratio of incorporated atoms over time—nanoparticle simulation The black line shows the ratio between all atoms in building blocks that are incorporated into a cluster compared to the number of atoms in the initial 40 building blocks (16 LCPA, 24 NP). The development of the largest cluster is shown by the blue area. A building block is counted as a member of a cluster if it is closer than 2.8 Å to it. The simulation setup is shown in Fig. 9

While the system with 16 LCPAs is already significantly larger than the basic setup from Fig. 7, it suffers from some limitations related to the limited box size and high concentration of the components. After the simulations, the final clusters connect to themselves over the box boundary surfaces (see Fig. 9b). These interactions could distort our results due to an additional constraint that will be entailed on the clusters. Therefore, in the second extended system we expanded the edges of the simulation box to 20 ns (see Fig. 11a). With more than 15 times the available volume, we can now simulate 120 LCPAs and 360 \({{\hbox {SiO}}_2}\) particles together. In this composition, a charge difference of 255 e was canceled out by chloride ions. For 47 ns 771,372 atoms were simulated. In Additional file 4, the simulation is displayed color-coded. Depending on the last cluster a building block belonged to, it is assigned a color. This assignment makes it easier to follow the development of the final clusters.

Fig. 11

Aggregation behavior of a larger polyamine-silica nanoparticle ensemble. a Snapshot of the initial configuration of 360 \({{\hbox {SiO}}_2}\) particles (gray; \(d\approx 10 - {20}\) Å) with 120 LCPA molecules (green) in a 20 nm cube of water. b Snapshot of the final stage of the simulation after 47 ns. The supercell used is highlighted with a white square surrounded by repeating neighboring supercells to visualize the periodic boundary conditions used in the MD simulations

Figure 12 displays the evolution of the generated clusters for the extended setup. In the beginning, free building blocks are incorporated fast into the clusters. After 5 ns just 5 % of all atoms are still available in the free components. Around the same time, the nearly linear growth of the five most significant clusters ends. After the initial phase that was driven by the incorporation of free building blocks, the growth is now bound to the interaction of the clusters themselves. This change can also be observed by the now sharp changes in the size of the tracked top five clusters in Fig. 12. As stated earlier, the combination of multiple clusters into a new larger one is hindered by the slower diffusion and the more limited number of adequate anchoring points. Despite these challenges, after 20 ns around 50 % of the available atoms are already incorporated into the five largest clusters. By the end of the simulation, this number grows to 78 %, demonstrating that the clusters are further stabilized. Figure 11b also demonstrates the fact that nearly all building blocks are included in a cluster at this point, as we find the areas between the now formed clusters empty. The shape of the main final clusters is approximately cylindrical with a length of 8 nm to 18 nm and a diameter between 2 nm and 4 nm.

Fig. 12

Ratio of incorporated atoms over time—extended nanoparticle simulation The black line shows the ratio between all atoms in the building blocks that are incorporated into a cluster compared to the number of atoms in the initial building blocks consisting of 480 components (120 LCPA, 360 NP). The five biggest clusters are stacked with the largest cluster on top in yellow and the fifth largest cluster at the bottom of the stack in violet. A building block is counted as a member of a cluster if it is closer than 2.8 Å to it. The simulation setup is shown in Fig. 11

Even this large simulation box introduces some artifacts that we would not expect in a more realistic scenario. In our simulations, the cluster starts to form a network-like structure. While thin water layers still separated the clusters, they clearly influence each other. Due to the nature of the periodic boundary conditions, these interactions equally distribute over the simulation space. Without this condition, the formation of a single spherical cluster would be expected based on the random attachments of new components. In a lower concentration setup, where the precursors would be replenished continuously, the linear growth of the clusters, like we see at the start of the simulations, would be more dominant. Also, the cluster–cluster interactions would play a less significant role.


We have carried out all-atom Molecular Dynamics simulations of the aggregation behavior of polyamines with phosphate ions and with various types of silica precursors. Concerning the aggregation of LCPAs in solution in the presence of phosphates, our results show that at a critical phosphate concentration of [Pi]/[a.g.]\(^+= 1.0\) the LCPA clusters display a tendency to remain structurally stable in an aqueous environment. Below this concentration, the solvation of phosphate ions by the water molecules dominates the dynamics, so that the aggregate fully dissolves. This behavior as a function of the Pi concentration points at the increasing role of electrostatic interactions with increasing Pi concentration, and it is in qualitative agreement with experimental results, see e.g., [2, 65]. However, a too large ionic concentration may jeopardize the structural stability of the aggregates due to an increase in the electrostatic repulsion forces. Moreover, our simulations also reveal, besides the Coulomb interactions, the existence of a hydrogen bond network across the phosphate ions, which also supports the conclusions of Ref. [38]. In the presence of polysilicic acid, the behavior of the aggregates upon immersion in aqueous solution is different for electrically neutral (dimers, tetramers) and negatively charged (octamers) species. For the former, the polyamine-silica aggregates rapidly dissolve in water after only a few nanoseconds simulation time. On the contrary, the aggregates built with the octamers mostly kept their structural integrity over the whole simulation time. This result shows again the dominant role played by electrostatic interactions in controlling the structural stability of the system. Hydrogen bond and van-der-Waals forces, dominant in the case of dimers and tetramers, do not seem strong enough to avoid the dissolution of the aggregates in water. To move one step further, we addressed in the second part of our study the aggregation of larger negatively charged silica nanoparticles mediated by LCPAs. Here, a clear tendency was found, showing the merging of many small clusters into only a few larger ones during the simulation time. These larger aggregates remain structurally stable over the time windows addressed by our simulations. In this case, the aqueous environment does not play any significant role in influencing the stability of the aggregates, which may be related to the full dominance of electrostatic interactions. These results are also in qualitative agreement with in vitro experiments [2] and highlight the potentially catalytic and stabilizing effect of LCPAs in biosilica aggregation.

Our simulations provide atomistic insight into the potentially dominant interactions in the early stages of biosilica formation. One major open problem, which cannot be addressed by classical molecular dynamics simulation approaches, is the very initial stage of silicic acid polycondensation, an issue related to the hypothesized role of polyamines in the formation of silica dimers via a proton exchange mechanism [8]. This effect, as well as bond formation, may play a role at the level of LCPA-dimer and LCPA-tetramer interactions; their treatment requires, however, a quantum chemical approach that goes beyond the scope of the present work.

Finally, we remark that a genuinely multi-scale approach to biosilica morphogenesis is still at its infancy. It will require the development of computational strategies to identify and transfer relevant physico-chemical information from the atomistic scale, as e.g., presented in the current study, to methodologies based either on strongly coarse-grained molecular dynamics or continuum models able to describe the time evolution and formation of biosilica patterns in diatoms or the more general context of biomineralization. In this respect, novel approaches involving machine learning techniques [71, 72] may prove to be of invaluable help in optimizing the related computational effort.

Availability of data and materials

Simulation data can be obtained upon request by the authors.



dihydrogen phosphate ion (\({{{\hbox {H}}}_2{\hbox {PO}}_4^-}\))


amino group


repeating units


disilicic acid (\({{{\hbox {Si}}}_2{{\hbox {O}}}_7{{\hbox {H}}}_6}\))


cyclo tetrasilicic acid (\({\hbox {Si}}_4{\hbox {O}}_{12}{\hbox {H}}_{8}\))


octahydroxyoctasilsequioxane-ion (\({{\hbox {Si}}_8{\hbox {O}}_{20}{\hbox {H}}_7^-}\))


long-chain polyamine


Van der Waals interaction


hydrogen bond interaction




molecular dynamics


  1. 1.

    Gordon R, Losic D, Tiffany MA, Nagy SS, Sterrenburg FAS. The Glass Menagerie: diatoms for novel applications in nanotechnology. Trends Biotechnol. 2009;27(2):116–27.

    CAS  Article  Google Scholar 

  2. 2.

    Sumper M, Brunner E. Learning from diatoms: nature’s tools for the production of nanostructured silica. Adv Funct Mater. 2006;16(1):17–26.

    CAS  Article  Google Scholar 

  3. 3.

    Kröger N, Poulsen N. Diatoms-from cell wall biogenesis to nanotechnology. Ann Rev Genet. 2008;42(1):83–107.

    CAS  Article  Google Scholar 

  4. 4.

    Poulsen N, Kröger N. Silica morphogenesis by alternative processing of silaffins in the diatom Thalassiosira pseudonana. J Biol Chem. 2004;279(41):42993–42999.

  5. 5.

    Brunner E, Richthammer P, Ehrlich H, Paasch S, Simon P, Ueberlein S, et al. Chitin-based organic networks: an integral part of cell wall biosilica in the diatom Thalassiosira pseudonana. Angew Chemie Int Ed. 2009;48(51):9724–7.

    CAS  Article  Google Scholar 

  6. 6.

    Sumper M, Brunner E, Lehmann G. Biomineralization in diatoms: characterization of novel polyamines associated with silica. FEBS Lett. 2005;579(17):3765–9.

    CAS  Article  Google Scholar 

  7. 7.

    Richthammer P, Börmel M, Brunner E, van Pée KH. Biomineralization in diatoms: the role of silacidins. ChemBioChem. 2011;12(9):1362–6.

    CAS  Article  Google Scholar 

  8. 8.

    Kröger N, Deutzmann R, Bergsdorf C, Sumper M. Species-specific polyamines from diatoms control silica morphology. Proc Natl Acad Sci. 2000;97(26):14133–14138.

  9. 9.

    Wenzl S, Hett R, Richthammer P, Sumper M. Silacidins: Highly acidic phosphopeptides from diatom shells assist in silica precipitation in vitro. Angew Chem Int Ed. 2008;47(9):1729–32.

    CAS  Article  Google Scholar 

  10. 10.

    Kröger N, Deutzmann R, Sumper M. Polycationic peptides from diatom biosilica that direct silica nanosphere formation. Science. 1999;286(5442):1129–1132.

  11. 11.

    Lechner CC, Becker CFW. Exploring the effect of native and artificial peptide modifications on silaffin induced silica precipitation. Chem Sci. 2012;3:3500–4.

    CAS  Article  Google Scholar 

  12. 12.

    Poulsen N, Sumper M, Kröger N. Biosilica formation in diatoms: characterization of native silaffin-2 and its role in silica morphogenesis. Proc Natl Acad Sci. 2003;100(21):12075–12080.

  13. 13.

    Bernecker A, Wieneke R, Riedel R, Seibt M, Geyer A, Steinem C. Tailored synthetic polyamines for controlled biomimetic silica formation. J Am Chem Soc. 2010;132(3):1023–31.

    CAS  Article  Google Scholar 

  14. 14.

    Wieneke R, Bernecker A, Riedel R, Sumper M, Steinem C, Geyer A. Silica precipitation with synthetic silaffin peptides. Org Biomol Chem. 2011;9:5482–6.

    CAS  Article  Google Scholar 

  15. 15.

    Jantschke A, Koers E, Mance D, Weingarth M, Brunner E, Baldus M. Insight into the supramolecular architecture of intact diatom biosilica from DNP-supported solid-state NMR spectroscopy. Angew Chem Int Ed. 2015;54(50):15069–73.

    CAS  Article  Google Scholar 

  16. 16.

    Kröger N, Lorenz S, Brunner E, Sumper M. Self-assembly of highly phosphorylated silaffins and their function in biosilica morphogenesis. Science. 2002;298(5593):584–586.

  17. 17.

    Perry CC, Keeling-Tucker T. Biosilicification: the role of the organic matrix in structure control. 2000;5:537–550.

  18. 18.

    Peguiron A, Colombi Ciacchi L, De Vita A, Kermode JR, Moras G. Accuracy of buffered-force QM/MM simulations of silica. J Chem Phys. 2015;142(6):064116.

    CAS  Article  Google Scholar 

  19. 19.

    Trinh TT, Jansen APJ, van Santen RA. Mechanism of oligomerization reactions of silica. J Phys Chem B. 2006;110(46):23099–106.

    CAS  Article  Google Scholar 

  20. 20.

    Zhang XQ, Trinh TT, van Santen RA, Jansen APJ. Mechanism of the initial stage of silicate oligomerization. J Am Chem Soc. 2011;133(17):6613–25.

    CAS  Article  Google Scholar 

  21. 21.

    Lenoci L, Camp PJ. Self-assembly of peptide scaffolds in biosilica formation: computer simulations of a coarse-grained model. J Am Chem Soc. 2006;128(31):10111–7.

    CAS  Article  Google Scholar 

  22. 22.

    Delle Piane M, Potthoff S, Brinker CJ, Colombi Ciacchi L. Molecular dynamics simulations of the silica-cell membrane interaction: insights on biomineralization and nanotoxicity. J Phys Chem C. 2018;122(37):21330–43.

    CAS  Article  Google Scholar 

  23. 23.

    Pérez-Sánchez G, Chien SC, Gomes JRB, Cordeiro MNDS, Auerbach SM, Monson PA, et al. Multiscale model for the templated synthesis of mesoporous silica: the essential role of silica oligomers. Chem Mater. 2016;28(8):2715–27.

    CAS  Article  Google Scholar 

  24. 24.

    Centi A, Jorge M. Molecular simulation study of the early stages of formation of bioinspired mesoporous silica materials. Langmuir. 2016;32(28):7228–40.

    CAS  Article  Google Scholar 

  25. 25.

    Jorge M, Milne AW, Sobek ON, Centi A, Pérez-Sánchez G, Gomes JRB. Modelling the self-assembly of silica-based mesoporous materials. Mol Simul. 2018;44(6):435–52.

    CAS  Article  Google Scholar 

  26. 26.

    Centi A, Manning JRH, Srivastava V, van Meurs S, Patwardhan SV, Jorge M. The role of charge-matching in nanoporous materials formation. Mater Horiz. 2019;15:1027–1033.!divAbstract.

  27. 27.

    Lenoci L, Camp PJ. Diatom structures templated by phase-separated fluids. Langmuir. 2008;24(1):217–23.

    CAS  Article  Google Scholar 

  28. 28.

    Schoeppler V, Gránásy L, Reich E, Poulsen N, de Kloe R, Cook P, et al. Biomineralization as a paradigm of directional solidification: a physical model for molluscan shell ultrastructural morphogenesis. Adv Mater. 2018;30(45):1803855.

    CAS  Article  Google Scholar 

  29. 29.

    Schoeppler V, Lemanis R, Reich E, Pusztai T, Gránásy L, Zlotnikov I. Crystal growth kinetics as an architectural constraint on the evolution of molluscan shells. Proc Natl Acad Sci. 2019;116(41):20388–20397.

  30. 30.

    Parkinson J, Brechet Y, Gordon R. Centric diatom morphogenesis: a model based on a DLA algorithm investigating the potential role of microtubules. Biochimica et Biophysica Acta (BBA)—Molecular Cell Research. 1999;1452(1):89–102.

    CAS  Article  Google Scholar 

  31. 31.

    Sumper M. A phase separation model for the nanopatterning of diatom biosilica. Science. 2002;295(5564):2430–2433.

  32. 32.

    Pohnert G. Biomineralization in diatoms mediated through peptide- and polyamine-assisted condensation of silica. Angew Chem Int Edn. 2002;41(17):3167–9.

    CAS  Article  Google Scholar 

  33. 33.

    Hildebrand M, Lerch SJL. Diatom silica biomineralization: parallel development of approaches and understanding. Semin Cell Dev Biol. 2015;46:27–35 (Biomineralisation and motorisation of pathogens).

    CAS  Article  Google Scholar 

  34. 34.

    Foo CWP, Huang J, Kaplan DL. Lessons from seashells: silica mineralization via protein templating. Trends Biotechnol. 2004;22(11):577–85.

    CAS  Article  Google Scholar 

  35. 35.

    Sumper M, Lorenz S, Brunner E. Biomimetic control of size in the polyamine-directed formation of silica nanospheres. Angew Chem. 2003;115(42):5350–3.

    Article  Google Scholar 

  36. 36.

    Sumper M. Biomimetic patterning of silica by long-chain polyamines. Angew Chem Int Edn. 2004;43(17):2251–4.

    CAS  Article  Google Scholar 

  37. 37.

    Brunner E, Lutz K, Sumper M. Biomimetic synthesis of silica nanospheres depends on the aggregation and phase separation of polyamines in aqueous solution. Phys Chem Chem Phys. 2004;6:854–7.

    CAS  Article  Google Scholar 

  38. 38.

    Lutz K, Gröger C, Sumper M, Brunner E. Biomimetic silica formation: analysis of the phosphate-induced self-assembly of polyamines. Phys Chem Chem Phys. 2005;7:2812–5.

    CAS  Article  Google Scholar 

  39. 39.

    Bridoux MC, Ingalls AE. Structural identification of long-chain polyamines associated with diatom biosilica in a Southern Ocean sediment core. Geochimica et Cosmochimica Acta. 2010;74(14):4044–57.

    CAS  Article  Google Scholar 

  40. 40.

    Abacilar M, Daus F, Haas C, Brückner SI, Brunner E, Armin Geyer A. Synthesis and NMR analysis of 13C and 15N-labeled long-chain polyamines (LCPAs). RSC Adv. 2016;6:93343–8.

    CAS  Article  Google Scholar 

  41. 41.

    Kröger N, Deutzmann R, Bergsdorf C, Sumper M. Species-specific polyamines from diatoms control silica morphology. PNAS. 2000;26(97):14133–14138. Available from:

  42. 42.

    Sumper M, Lorenz S, Brunner E. Biomimetic control of size in the polyamine-directed formation of silica nanospheres. Angew Chem Int Ed Engl. 2003;42(42):5192–5.

    CAS  Article  Google Scholar 

  43. 43.

    Sumper M. Biomimetic patterning of silica by long-chain polyamines. Angew Chem Int Ed Engl. 2004;17(43):14133–8.

    CAS  Article  Google Scholar 

  44. 44.

    Lutz K, Gröger C, Manfred Sumper M, Brunner E. Biomimetic silica formation: analysis of the phosphate-induced self-assembly of polyamines. Phys Chem Chem Phys. 2005;7(14):2812–2815.!divAbstract.

  45. 45.

    Campbell NA, Reece JB. Biology (7h ed.) Benjamin Cummings, San Francisco, p 65 ISBN 0-8053-7171-0. 2005.

  46. 46.

    Belton DJ, Deschaume O, Perry CC. An overview of the fundamentals of the chemistry of silica with relevance to biosilicification and technological advances. FEBS J. 2012;279:1710–20.

    CAS  Article  Google Scholar 

  47. 47.

    Cruz-Chu ER, Aksimentiev A, Schulten K. Water-silica force field for simulating nanodevices. J Phys Chem B. 2006;110(43):21497–508.

    CAS  Article  Google Scholar 

  48. 48.

    Huff NT, Demiralp E, Çagin T, Goddard WA III. Factors affecting molecular dynamics simulated vitreous silica structures. J Non Cryst Solids. 1999;253(1–3):133–42.

    CAS  Article  Google Scholar 

  49. 49.

    Emami FS, Puddu V, Berry RJ, Varshney V, Patwardhan SV, Perry CC, et al. Force field and a surface model database for silica to simulate interfacial properties in atomic resolution. Chem Mater. 2014;26(8):2647–58.

    CAS  Article  Google Scholar 

  50. 50.

    Phillips JC, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, et al. Scalable molecular dynamics with NAMD. J Comput Chem. 2005;26:1781–802.

    CAS  Article  Google Scholar 

  51. 51.

    Hornak V, Abel R, Okur A, Strockbine B, Roitberg A, Simmerling C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins. 2006;3(65):712–25.

    Article  Google Scholar 

  52. 52.

    Heinz H, Lin TJ, Mishra RK, Emami FS. Thermodynamically consistent force fields for the assembly of inorganic, organic, and biological nanostructures: the INTERFACE force field. Langmuir. 2013;6(29):1754–65.

    CAS  Article  Google Scholar 

  53. 53.

    Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, et al. General atomic and molecular electronic structure system. J Comput Chem. 1993;14:1781–1802.

  54. 54.

    Vanquelef E, Simon S, Marquant G, Garcia E, Klimerak G, Delepine J, et al. Scalable molecular dynamics with NAMD. Nucleic Acids Res. 2011;39:1781–802.

    Article  Google Scholar 

  55. 55.

    Dupradeau FY, Pigache A, Zaffran T, Savineau C, Lelong R, Grivel N, et al. The R.E.D. tools: advances in RESP and ESP charge derivation and force field library building. Phys Chem Chem Phys. 2010;12:7821–39.

    CAS  Article  Google Scholar 

  56. 56.

    Ryckaert JP, Ciccotti G, Berendsen HJC. Numerical-integration of cartesian equations of motion of a system with constraints—molecular-dynamics of N-alkanes. J Comput Phys. 1977;3(23):327–41.

    Article  Google Scholar 

  57. 57.

    Darden T, York D, Pedersen L. Particle mesh Ewald: An N-log(N) method for Ewald sums in large systems. J Chem Phys. 1993;98:10089–92.

    CAS  Article  Google Scholar 

  58. 58.

    Haas KR, Yang H, Chu JW. Fisher information metric for the Langevin equation and least informative models of continuous stochastic dynamics. J Chem Phys. 2013;139:121931.

    Article  Google Scholar 

  59. 59.

    Humphrey W, Dalke A, Schulten K. VMD: Visual molecular dynamics. J Mol Graph. 1996;14(1):33–8.

    CAS  Article  Google Scholar 

  60. 60.

    Stone J. An efficient library for parallel ray tracing and animation. Rolla: Computer Science Department, University of Missouri-Rolla; 1998.

    Google Scholar 

  61. 61.

    Michaud-Agrawal N, Denning EJ, Woolf TB, Beckstein O. MDAnalysis: a toolkit for the analysis of molecular dynamics simulations. J Comput Chem. 2011;32(10):2319–27.

    CAS  Article  Google Scholar 

  62. 62.

    Gowers R, Linke M, Barnoud J, Reddy T, Melo M, Seyler S, et al. MDAnalysis: a python package for the rapid analysis of molecular dynamics simulations. In: Proceedings of the 15th python in science conference. Scipy; 2016. p. 98–105. Available from:

  63. 63.

    Eckert H. Molani: Automated molecular animator; 2020.

  64. 64.

    Brunner E, Lutz K, Sumper M. Biomimetic synthesis of silica nanospheres depends on the aggregation and phase separation of polyamines in aqueous solution. Phys Chem Chem Phys. 2004;6:854–857. Available from:!divAbstract.

  65. 65.

    Agazzi LM, Herrera SE, Cortez ML, Marmisolle WA, von Bilderling C, Pietrasanta LI, et al. Continuous assembly of supramolecular polyamine–phosphate networks on surfaces: preparation and permeability properties of nanofilms. Soft Matter. 2019;15:1640–1650. Available from:!divAbstract.

  66. 66.

    Sumper M, Brunner E. Learning from diatoms: nature’s tools for the production of nanostructured silica. Adv Funct Mater. 2006;16(1):17–26.

    CAS  Article  Google Scholar 

  67. 67.

    Schmid AM, Volcani BE. Wall morphogenesis in coscinodiscus wailesii gran and angst. I. Valve morphology and development of its architecture. J Phycol. 1983;19(4):387–402.

    Article  Google Scholar 

  68. 68.

    Di Luccia A, Picariello G, Iacomino G, Formisano A, Paduano L, D’Agostino L. The in vitro nuclear aggregates of polyamines. FEBS J. 2009;276(8):2324–35.

    CAS  Article  Google Scholar 

  69. 69.

    Savitzky A, Golay MJE. Smoothing and differentiation of data by simplified least squares procedures. Anal Chem. 1964;36(8):1627–39.

    CAS  Article  Google Scholar 

  70. 70.

    Berendsen HJC, van der Spoel Dv, Drunen R. GROMACS: a message-passing parallel molecular dynamics implementation. Comput Phys Commun. 1995;91:45–56.

    Article  Google Scholar 

  71. 71.

    Wang J, Olsson S, Wehmeyer C, Pérez A, Charron NE, de Fabritiis G, et al. Machine learning of coarse-grained molecular dynamics force fields. ACS Cent Sci. 2019;5(5):755–67.

    CAS  Article  Google Scholar 

  72. 72.

    Noé F, Tkatchenko A, Müller KR, Clementi C. Machine learning for molecular simulation. Annu Rev Phys Chem. 2020;71(1):361–90.

    CAS  Article  Google Scholar 

Download references


The authors thank the Deutsche Forschungsgemeinschaft (FOR 2038: Nanopatterned Organic Matrices in Biological Silica Mineralization) for financial support and the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for computational resources. The authors also thank Eike Brunner (TU Dresden, Germany), N. Kröger (BCube Dresden, Germany), Armin Geyer (Philipps University of Marburg, Germany), and Miguel Jorge (University of Strathclyde, UK) for very fruitful and interesting discussions.


Deutsche Forschungsgemeinschaft: Nanopatterned Organic Matrices in Biological Silica Mineralization Research Unit FOR 2038). The DFG only provided funding for the positions of H. E. and M. M. The DFG did not contribute to the design of the study and collection, analysis, and interpretation of data or in writing the manuscript. Open Access Funding was provided by the Publication Fund of the TU Dresden.

Author information




H.E. and M.M. carried out the computer simulations, M.B., R.G., A.D., and G.C. wrote and corrected the manuscript. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Rafael Gutierrez.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information


Additional file 1. Solving random dimer cluster in water. Simulation of a pre-aggregated cluster composed of ninefold protonated polyamines, negatively phosphate ions, and electrically neutral polysilicic acid dimers upon immersion in aqueous environment. The aggregate quickly destabilizes (within the first 10 ns) and finally fully dissolves in water.


Additional file 2. Solving random octamer cluster in water. Simulation of a pre-aggregated cluster composed of ninefold protonated polyamines, negatively charged phosphate ions, and negatively charged polysilicic acid octamers upon immersion in aqueous environment. Although during the first nanoseconds of the simulation most of the phosphate ions (red) and few octamers (gray) detach from the aggregate, it nevertheless largely preserves its structural integrity after 40 ns simulation time.


Additional file 3. Clustering of a single LCPA with two silica nanoparticles in water. Simulation of two SiO2 particles (d ≈ 17 Å) with one LCPA molecule in water for 20 ns. Already after 10 ns the system clusters due to the strong electrostatic interactions between the ninefold protonated moleculae and the negatively charged nanoparticles.


Additional file 4. Large scale aggregation of LCPA and silica nanoparticles. Simulation of 360 SiO2 particles (d ≈ 10–20 Å) with 120 LCPA molecule in water for 47 ns.

Additional file 1. Solving random dimer cluster in water. Simulation of a pre-aggregated cluster composed of ninefold protonated polyamines, negatively phosphate ions, and electrically neutral polysilicic acid dimers upon immersion in aqueous environment. The aggregate quickly destabilizes (within the first 10 ns) and finally fully dissolves in water.

Additional file 2. Solving random octamer cluster in water. Simulation of a pre-aggregated cluster composed of ninefold protonated polyamines, negatively charged phosphate ions, and negatively charged polysilicic acid octamers upon immersion in aqueous environment. Although during the first nanoseconds of the simulation most of the phosphate ions (red) and few octamers (gray) detach from the aggregate, it nevertheless largely preserves its structural integrity after 40 ns simulation time.

Additional file 3. Clustering of a single LCPA with two silica nanoparticles in water. Simulation of two SiO2 particles (d ≈ 17 Å) with one LCPA molecule in water for 20 ns. Already after 10 ns the system clusters due to the strong electrostatic interactions between the ninefold protonated moleculae and the negatively charged nanoparticles.

Additional file 4. Large scale aggregation of LCPA and silica nanoparticles. Simulation of 360 SiO2 particles (d ≈ 10–20 Å) with 120 LCPA molecule in water for 47 ns.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Eckert, H., Montagna, M., Dianat, A. et al. Exploring the organic–inorganic interface in biosilica: atomistic modeling of polyamine and silica precursors aggregation behavior. BMC Mat 2, 6 (2020).

Download citation


  • Diatoms
  • Atomistic molecular dynamics
  • Biosilica formation
  • Long-chain polyamines
  • Aggregation