 Research article
 Open Access
 Published:
Continuum modelling of structure formation of biosilica patterns in diatoms
BMC Materials volume 2, Article number: 12 (2020)
Abstract
Formation of regularly structured silica valves of various diatom species is a particularly fascinating phenomenon in biomineralization. Intensive investigations have been devoted to elucidate the formation mechanisms of diatom valve structures. Phaseseparation of speciesspecific organic molecules has been proposed to be involved in pattern formation, where the evolving organic molecule structures serve as template for silica formation. In the present work, using a continuum approach, we investigate the conditions under which silica structures of high regularity can develop within a phase separation model. In relation to previously reported in vitro experiments of silica formation, which revealed the important role of phosphate ions in the selfassembly of organic molecules, we propose a model where phase separation is coupled with a chemical reaction. We analyze the impact of the reaction of phosphate ions with organic molecules on the appearing morphology of the organic template. Two and threedimensional simulations of the development of regular stationary patterns are presented. The influence of a confined geometry and an interaction of organic molecules with the walls on pattern formation is also addressed. We expect that our approach will be relevant for experimental studies aiming at inducing structure formation under controlled in vitro conditions.
Introduction
The formation of various beautiful silica structures of diatoms is a very fascinating example of biomineralization. In recent years, considerable progress has been made in elucidating the chemical and biological processes, which are involved in the silica morphogenesis of diatom valves (cf. for example [1,2,3,4,5]). Summaries of very recent findings have been given in [6, 7]. So far, only relatively few works have dealt with the very challenging theoretical modeling of the morphogenesis of diatom valves [8,9,10,11,12,13]. Early models were based on diffusionlimited aggregation of silica nanoparticles, which are imported into the silica deposition vesicle (SDV) at its leading edge [9, 11]. Corresponding computer simulations were able to describe the formation of the radial rib structure of valves. On the basis of diffusionreaction mechanisms, the formation of various nanosized pore patterns within larger pores of silica valves has also been addressed in [12, 13]. A special computer evolution model of the morphogenesis of raphid pennate diatoms was presented in [8]. Another modeling approach was based on the hypothesis that silica formation is guided by a template, which is formed by the selfassembly of organic components [1, 3, 4, 10, 14]. In this context, phase separation of organic molecules within the SDV was proposed as one formation mechanism of the patterned templates [14]. In a corresponding theoretical study [10], numerous twodimensional patterns were obtained, which showed great similarity with observed valve structures. In those computer simulations, phase separation was influenced by an additional local field arising from the preexisting silica costae within the base layer. Depending on the choice of this regular prepatterning field, different highly symmetric valve structures could be achieved.
Extensive experimental investigations have revealed that speciesspecific organic molecules are involved in the formation of the biosilica structures in diatoms, for example long chain polyamines, silaffins and silacidins [4, 7]. Due to the large complexity of in vivo studies, numerous in vitro experiments of the silification of aggregates of organic molecules have been performed to elucidate the mechanisms of biosilica formation [2, 15, 16]. It was found that the amount, size, and form of silica precipitates sensitively depend on the molecular structure of the synthesized polyamines. Further, the influence of phosphate ions on the precipitation was analyzed. In particular, the analysis in [2] showed that silica precipitation occurred beyond a certain threshold pH which depends on the phosphate concentration. With increasing phosphate concentration, the diameter of spherical aggregates of the organic molecules also increased. Those aggregates served as template for the silica formation after adding silicic acid molecules to the solution. The diameter of the arising silica spheres was about 100 to 600 nm [2]. Interestingly, the size of the polyamine aggregates could be controlled also at fixed polyamine/phosphate ratio by changing the pH value. In [15], besides hollow spheres of about 20 to 200 nm size, filamentlike silica structures were also observed. The described silification of spherical organic templates does not explain, however, the formation of the various pore structures of diatom valves. To obtain a porous silica structure, an appropriate mechanism was proposed in [2]. Instead of monosilicic acid as silicon source, a sol of stabilized silica nanoparticles was used for the formation of silica structures. The nanoparticles arranged around the organic aggregates, so that a networklike silica morphology was formed. The longterm goal of in vitro approaches would be the creation of artificial SDVs, where biosilica formation can be controlled by varying key parameters such as the pH and the concentration of the different participating components (organic and inorganic). Achieving this goal would be not only of great relevance for understanding the underlying mechanisms leading to biosilicification in vivo, but it would also trigger farreaching nanotechnological applications of diatoms [17].
Inspired by the previously described in vitro experiments and by the proposals in [2, 14], we will investigate in this work phase separation of organic molecules as one possible mechanism to form organic templates for the formation of silica structures. To avoid the permanent coarsening of phase regions in conventional phase separation, we propose here a model where the phase separation is coupled with a chemical reaction [18], an approach which can be justified on the basis of the in vitro experiments of silica precipitation. In particular, we focus on the role of phosphate anions in the selfassembly of longchain polyamines (LCPAs). Besides the experimental studies mentioned above, the importance of phosphate ions is further highlighted by experiments on silaffins, which showed that silaffin 1A (in contrast to natSil1A) was not able to promote silica precipitation in the absence of phosphates [19]. This was traced back to the absence in silaffin 1A of phosphorylation of the hydroxyamino acids of the peptide backbone. The choice of LCPAs as a reference organic component is motivated by their demonstrated influence on silica aggregation in vitro [20,21,22,23] as well as by their broad dispersion into the silica phase, shown recently in [24]. The latter fact suggests their relevance for the formation of the organic–inorganic interface under in vivo conditions.
In view of these findings, we suggest that the reaction of organic molecules with phosphate ions has an important impact on the phase separation of the organic component and, therefore, it should be taken into account in the theoretical modeling of this process. Our aim was to find a mechanism how regular stationary phase patterns develop.
The paper is outlined as follows. First, we present our mathematical model of the phase separation of organic molecules coupled to a chemical reaction. Then, by means of a linear stability analysis of the model equations, we show under which conditions spinodal decomposition of a uniform state occurs. As a result of our numerical simulations, various stationary phase patterns of high regularity are obtained and the impact of different model parameters on the emerging patterns is analyzed.
Model
Our theoretical investigation of the aggregation of organic molecules is based on the following approach. To study the role of phosphate ions in the aggregation process, we consider the simple generic reaction:
where component A corresponds to organic molecules which are unable to phase separate. We hypothesize that by binding of phosphate ions (component B) to those organic molecules, the organic component C is formed,which exhibits a higher degree of phosphorylation so that phase separation of component C becomes possible due to its reduced charge compared to component A. More explicitly, the above reaction may adopt the following form :
To describe the phase separation process, we use the wellknown CahnHilliard model [10, 25, 26]. We start from the Gibbs energy functional:
where \(\phi (\vec {r})\) is the volume fraction of the phase separating component and \(g(\phi )\) is the Gibbs energy density of mixing. The gradient energy coefficient \(\kappa \) will be assumed to be constant [27]. Together with the constraint \({\frac{1}{V}\int {\rm{d}}V \phi (\vec {r}) = \bar{\phi }}\) and the continuity equation
with \(M\) being the mobility and \(\mu =\delta G/\delta \phi \) the chemical potential, the evolution equation for the phase field \(\phi (\vec {r},t)\) is obtained as [10, 26]:
For simplicity, the mobility was approximated as a constant. For the Gibbs energy density, we adopted a generic asymmetric double well potential according to the FloryHuggins model:
where \(g_0\) is a prefactor, \(N\) is the degree of polymerization, and \(\chi \) is the FloryHuggins interaction parameter [27]. Spinodal decomposition occurs for volume fractions where \(g''(\phi )<0\). In the following calculations, we used the parameters \(N=10\) and \(\chi =1.3\). For this example, spinodal decomposition takes place for volume fractions \(0.066<\phi <0.588\). Introducing the length unit \(l=\sqrt{\kappa /g_0}\) and time unit \(\tau =l^2/(Mg_0)\), Eq. 4 can be written in dimensionless form as:
where \(\hat{t}=t/\tau \), \(\hat{\nabla }=l\nabla \), and \(\hat{g}=g/g_0\). The unit length is of the order of magnitude of the interface width between the separated phases (cf. for example Fig. 3 below). The evolution of the concentrations of components A and B is described by the following reaction–diffusion equations [18]:
where \(c_{\rm{A}}\), \(c_{\rm{B}}\), and \(c_{\rm{C}}\) are the concentrations (with dimension \(\hbox {m}^{3}\)) of the corresponding components. The reaction velocity of components A and B is described by the coefficient \(\alpha \) and the backward reaction (dissociation of component C) by the coefficient \(\beta \). The diffusion coefficients \(D_{\rm{A}}\) and \(D_{\rm{B}}\) generally depend on the concentrations \(c_{\rm{A}}\) and \(c_{\rm{B}}\), on the volume fraction \(\phi \) [28, 29], and on the charge state of the polymer [30]. In fact, variations over several orders of magnitude can be found. Due to the lack of specific information on the studied polyaminephosphate ionwater system, we will assume in the following that these kinetic coefficients are constants in a first analysis. A rough orderofmagnitude estimation of the reactions rates \(\alpha \) and \(\beta \) is provided in the Supplementary Material. Since we decided to describe phase separation of component C in Eq. 6 in terms of the volume fraction of component C, Eqs. 7 and 8 have to be properly adapted. The volume fraction of C is related to the concentration by \(\phi =\Omega _{\rm{C}}c_{\rm{C}}\) with \(\Omega _{\rm{C}}\) as the molecular volume of component C. Thus, Eqs. 7 and 8 are multiplied by \(\Omega _{\rm{C}}\) and by using the above length and time units, one obtains:
with \(\hat{c}_{\rm{A}} = \Omega _{\rm{C}}c_{\rm{A}}\), \(\hat{c}_{\rm{B}} = \Omega _{\rm{C}}c_{\rm{B}}\), \(\hat{D}_{\rm{A}} = D_{\rm{A}}\tau /l^2\), \(\hat{D}_{\rm{B}} = D_{\rm{B}}\tau /l^2\), \(\hat{\alpha }=\alpha \tau /\Omega _{\rm{C}}\), and \(\hat{\beta }=\beta \tau \). Note that concentrations \(c_{\rm{A}}\) and \(c_{\rm{B}}\) are scaled by the molecular volume of component C. Adding the reaction terms in Eq. 6 finally yields:
For the sake of brevity, we will omit the roof on the scaled variables in the following.
It should be noted that the chosen expression for the Gibbs energy in Eq. 5 is valid for a polymersolvent system. Actually, an interaction of component C with components A and B should be also taken into account (cf. for example [27]). However, to simplify the situation in this first analysis, we will omit this interaction. In the limit of a small dissociation constant \(\beta \), the volume fractions of components A and B will be small compared to that of component C, in particular to that of component B, since the phosphate ions have a much smaller molecular volume than the organic component.
We further consider here the case that the reaction and diffusion processes are fast compared to the phase separation. Then, starting from an initial condition \(c_{\rm{A}} = c_{\rm{A,0}}\), \(c_{\rm{B}} = c_{\rm{B,0}}\) and \(\phi =0\), a reaction equilibrium \(\alpha c_{\rm{A,e}}c_{\rm{B,e}}=\beta \phi _{\rm{e}}\) will be quickly established. Together with the particle conservation laws \(c_{\rm{A,0}} = c_{\rm{A,e}} + \phi _{\rm{e}}\) and \(c_{\rm{B,0}} = c_{\rm{B,e}} + \phi _{\rm{e}}\), the equilibrium concentration of component C is obtained as:
Phase separation will not occur if this equilibrium concentration is outside the spinodal region of the Gibbs energy (\(g''(\phi )<0\)). The system of nonlinear partial differential equations 9–11, which describe the coupled evolution of the three components, has been solved numerically. In case of periodic boundary conditions, we mainly used an efficient semiimplicit Fourierspectral method [31], which we implemented in the software MATLAB R2019b. For comparison and for simulations within a confined geometry, we applied also the finite element software Comsol Multiphysics 5.2 [32]. Our model of the phase separation is characterized by the four model parameters \(D_{\rm{A}}\), \(D_{\rm{B}}\), \(\alpha \) and \(\beta \), as well as by the initial concentrations \(c_{\rm{A,0}}\) and \(c_{\rm{B,0}}\). In view of the great number of parameters, we analyze in a first step, under which conditions a uniform system becomes unstable and spinodal decomposition can occur.
Stability analysis
We start from a spatially uniform system with all concentrations at their equilibrium values, i.e. \(\alpha c_{\rm{A,e}}c_{\rm{B,e}}=\beta \phi _{\rm{e}}\), and analyze whether small fluctuations will unstably grow or decay. Inserting \(c_{\rm{A}} = c_{\rm{A,e}} + \delta c_{\rm{A}} \), \(c_{\rm{B}} = c_{\rm{B,e}} + \delta c_{\rm{B}} \), and \(\phi = \phi _{\rm{e}} +\delta \phi \) into Eqs. 9 to 11 and expanding up to linear order in the fluctuations, we obtain:
Inserting the ansatz \(\delta c_{\rm{A}} = a_{\rm{A}}\exp {(\omega t)}\exp {(iqx)}\), \(\delta c_{\rm{B}} = a_{\rm{B}}\exp {(\omega t)}\exp {(iqx)}\), \(\delta \phi = a_{\phi }\exp {(\omega t)}\exp {(iqx)}\) into Eqs. 13–15 yields the following equations for the amplitudes which, for brevity, are represented in matrix form:
This linear system of homogeneous equations has a solution for vanishing of the determinant of the coefficient matrix. This condition leads to a cubic equation for the frequency as a function of the wave vector. In the considered parameter region, the three solutions of the equation are real. Two of them are negative at all wave vectors. One solution, which exhibits also positive values, is referred to as the dispersion relation \(\omega (q)\). Spinodal decomposition of component C will occur for growing fluctuations, i.e. for qvalues for which \(\omega (q)>0\). In Fig. 1, examples of the dispersion relation are shown for different dissociation constants \(\beta \). In general, the frequency becomes negative above a critical wave vector \(q_0\) with \(\omega (q_0)=0\). Thus, fluctuations with wavelengths shorter than \(\lambda _{0}=2\pi /q_0\) will decay in time. The wavelength of the fastest growing Fourier mode is given by the maximum of \(\omega (q)\) at the wave vector \(q_{\text {max}}\). The corresponding wavelength \(\lambda _{max}=2\pi /q_{\text {max}}\) determines the characteristic size of the initially emerging phase regions.
In the limit of vanishing dissociation (\(\beta =0\)), the frequency \(\omega (q)\) is positive for all wave vectors \(q<q_0\). Thus, also fluctuations with very long–wavelengths will grow. As a consequence, phase regions will permanently coarsen in the nonlinear growth regime, which hinders the desired formation of regular stationary phase patterns. Interestingly, with increasing dissociation constant \(\beta \), the frequency within a certain region at small wave vectors becomes negative. This means that the corresponding longwavelength fluctuations are suppressed. As it will be shown below, the existence of such a suppressed wave vector region is of great importance for the establishment of regular stationary phase patterns, since it inhibits the permanent coarsening of phase regions. With further increase of the dissociation constant, the frequency \(\omega (q)\) becomes negative at all wave vectors and, hence, spinodal decomposition cannot occur. We remark that a similar dispersion relation would be obtained for the case where the phase separating component exhibits electrostatic repulsion. In this case, the Gibbs energy in Eq. 2 would need to be complemented by a term describing Coulomb interactions. The strength of this repulsion contribution has a similar effect on the dispersion relation as the dissociation constant in the present case.
In our stability analysis in Fig. 1a, b, we have given the initial concentrations \(c_{\rm{A,0}}\) and \(c_{\rm{B,0}}\) of components A and B (and \(\phi =0\)) and assumed a rapid establishment of reaction equilibrium according to Eq. 12. The change of the dispersion relation with increasing dissociation constant is shown in Fig. 1a. The growth velocity \(\omega _{\text {max}}\) diminishes with increasing dissociation constant, whereas the corresponding wavelength does not change. The impact of the initial concentration \(c_{\rm{B,0}}\) (phosphate ions) on the dispersion relation, while fixing the other parameters, is shown in Fig. 1b. The curves demonstrate that a sufficiently high phosphate concentration is necessary to enable spinodal decomposition by reaching a wave vector region with \(\omega (q)>0\). The wavelength of the fastest growing Fourier mode decreases slightly with increasing initial concentration \(c_{\rm{B,0}}\). The present stability analysis clearly characterizes only the very early stage of phase separation. The subsequent nonlinear concentration evolution has to be computed numerically. However, this analysis is of great benefit for finding out those parameter sets, where the formation of regular stationary phase patterns can be expected.
Numerical analysis of phase separation
1D simulations
To illustrate characteristic features of the nonlinear phase separation process, we first present results of numerical calculations in one dimension. As initial condition, the concentrations \(c_{\rm{A,0}}\) and \(c_{\rm{B,0}}\) of components A and B with small fluctuations of the order of 1 % are given. These initial fluctuations quickly decay due to diffusion processes. For component C, a very small initial concentration \(\phi =0.002\) is given. As a boundary condition, we chose periodic boundary conditions. Figure 2 shows an example of the concentration evolution at two times. At an intermediate time, several peaks appear. Subsequently, many peaks vanish and a state with two remaining peaks is achieved. Compared to Fig. 2b, the peaks further move up to a stationary state with a peak spacing of 40, i.e. half the cell size of \(\hbox {L} = 80\). Doubling the cell size to \(\hbox {L} = 160\) leads to a stationary pattern with five equidistant peaks with a spacing of 32 (cf. Supplementary Material). Further cell doubling results in nine peaks with a spacing of 35.56. Clearly, the periodicity length of the arising pattern is affected by the simulation cell size, which should be large compared to characteristic lengths of the pattern.
For demonstration purposes, comparatively small diffusion coefficients were chosen in the example in Fig. 2. It was further assumed that the diffusion coefficient of component B (phosphate ions) is considerably larger than that of component A (organic molecules). To ensure that the diffusion processes are fast compared to the phase separation process, in the following examples the diffusion coefficients \(D_{\rm{A}}={100}\) and \(D_{\rm{B}}={1000}\) for components A and B were chosen. (See the Supplementary Material for other simulations using smaller values of the diffusion coefficients. As a tendency, we found that for smaller diffusion coefficients the reaction and dissociation rates have also to be diminished to obtain regular patterns.) Furthermore, the reaction constant \(\alpha =10\) was kept fixed. The strong impact of the dissociation constant \(\beta \) on the concentration evolution is shown in Fig. 3. For a dissociation constant of \(\beta =0.15\), a regular arrangement of concentration peaks is quickly established, whereas without dissociation (\(\beta =0\)), a coarsening of phase regions occurred. The concentration profile for \(\beta =0.15\) in Fig. 3 is obviously stationary, since it did not change within our calculation time of up to \(t=10^{6}\). Note also the smaller difference between maximum and minimum concentration values in Fig. 3a, compared to the case without dissociation. The wavelength of the regular concentration profile for \(\beta =0.15\) in Fig. 3 roughly corresponds to the wave vector at the maximum of the dispersion relation in Fig. 1. The present onedimensional calculations reflect already important tendencies on how the choice of model parameters may affect the phase separation process in higher dimension.
2D simulations
In previous modelling studies on the structure development in diatoms, twodimensional (2D) models of pattern formation were considered. Besides computational simplicity, the restriction to two dimensions is related to the flat geometry of the silica deposition vesicle. Accordingly, we present in the following results of 2D simulations for the concentration evolution of component C. As boundary condition for the simulations in Fig. 4, we used periodic boundary conditions in all directions. Depending on the parameter choice, irregular patterns showing strong coarsening of phase regions (Fig. 4) are obtained, or very regular stationary patterns (Fig. 5). The regular pattern shows a hexagonal arrangement of aggregates of organic molecules of equal size. The orientation of this hexagonal array is to a certain degree accidental due to the random fluctuations of the initial concentration of components A and B. Further, as a finite size effect of our limited simulation cell, this pattern is influenced by the used periodic boundary conditions. By adding monosilicic acid to a solution with organic aggregates in in vitro experiments, silification of those aggregates was observed [15]. On the other hand, by adding a sol of silica nanoparticles, silica could form around the organic aggregates [2]. Corresponding to the regular aggregates in Fig. 5, a hexagonal pore structure would be obtained.
Figure 6 shows patterns of organic aggregates for higher initial concentrations of components A and B. The arising stationary patterns strongly depend on the value of the dissociation constant. Typically, a dense network of organicrich phase of component C with circular pores of organicpoor phase develops. The size of the pores increases with decreasing dissociation constant. Interestingly, at certain values of the dissociation constant, also lamellar (labyrinthlike) phase patterns were obtained. It is remarkable that those stripes are stationary despite their irregular structure.
Confinement effects. The phase patterns shown in Figs. 4, 5 and 6 were calculated by using periodic boundary conditions in order to minimize artificial finite size effects owing to the limited size of the simulation cell. In diatoms, the aggregation of organic material takes place within the SDV (the same may hold in the case of artificially built SDVs). This confinement is rather large compared to the characteristic lengths of the observed pore structures of valves of diatoms. Thus, it seems justified to consider only a small part of the SDV and to assume periodic boundary conditions in our above calculations. However, in the case of the observed hierarchical pore structure in diatom C. wailesii [14], preexisting large pores represent confinements for the formation of smaller pores. Pattern formation in preexisting pores was subject of investigations in [12, 13], which were based on a reaction–diffusion mechanism. We also performed simulations within a confined geometry based, however, on the phase separation mechanism. The patterns shown in Fig. 7 are computed within a circular simulation cell supposing zero flux boundary conditions. The obtained patterns are stationary within the simulation time of typically \(t=10^{5}\). Strictly, it cannot be excluded that some changes occur at very long time. Most of the patterns show a high regularity, depending on the initial concentration of components A and B and on the dissociation constant. If the size of the characteristic features of the patterns is small compared to the simulation cell, the patterns are only little influenced by the confinement. With decreasing dissociation constant, the feature size increases and the evolving patterns are strongly influenced by the circular confinement. Importantly, the formation of hexagonal arrays of small aggregates is observed within large pores as it was suggested by Sumper in [14]. At higher initial concentrations of components A and B, organicrich networks with circular pores or lamellar structures develop. Simulations with a \(\phi \)dependent mobility have also been carried out and are shown in Fig. S6 of the Supplementary Material.
Boundary interaction. A recent experimental study [33] focused on the influence of the boundary of the SDV on the aggregation of biomolecules. In particular, the interaction of longchain polyamines with lipid membranes was analyzed. Motivated by these investigations, we included in our modeling an attractive interaction of component A (organic molecules) with the boundary of the confinement. A corresponding simulation is shown in Fig. 8. To mimic the interaction with the boundary (left boundary in Fig. 8), we supposed an additional field \(h(x)=h_0\exp (x/\xi )\), acting on molecules A, and added a corresponding \(D_{\rm{A}} c_{\rm{A}} {\rm{d}}h(x)/{\rm{d}}x\) term on the right–hand side of the evolution equation for component A in Eq. 9. The field rapidly decays with increasing distance from the boundary. As a consequence of the enrichment of component A at the attracting boundary, a lamellar phase structure parallel to the boundary develops initially. Subsequently, the lamellae decompose into circular aggregates, forming a hexagonal array with defined orientation.
3D simulations
Microscopic images of the threedimensional valve structures show great differences between the proximal and distal sides of diatom valves. Thus, it is desirable to analyze the valve formation also by 3D computer simulations. In the following analysis, we have chosen a coordinate system with zdirection perpendicular to the flat SDV. Periodic boundary conditions were supposed only in x and ydirections, whereas at the boundaries in zdirection the normal diffusion flux for all components is required to vanish. Figure 9 shows examples for the establishment of stationary aggregates of component C. The aggregates are represented by surface plots of the concentration as well as by plots of the isoconcentration surface. The chosen kinetic parameters in Fig.9a, b are equal, only the thickness of the simulation cell differs slightly. For the example in Fig. 9a cylindrical aggregates developed, whereas in Fig. 9b halfsphere aggregates at the bottom and top face of the simulation cell were formed.
Conclusions
Following the hypothesis that the formation of silica structures of diatoms may be guided by an organic template, we studied the formation of aggregates of organic molecules by means of a phase separation model. This general hypothesis has received some additional support from extensive in vitro experiments, which have demonstrated the role of the organic components (in conjunction with the presence of multivalent anions) in catalyzing silica precipitation. In contrast to previous work [10], we did not assume an additional prepatterning field influencing the formation of highly regular patterns of the organic component. Instead, to avoid the permanent phase coarsening in common phase separation and to obtain patterns of high regularity, we coupled the phase separation process with a chemical reaction [18], where the phase separating component could dissociate into two components. Owing to this dissociation, regular stationary patterns of aggregated organic molecules could be obtained. We have demonstrated that the geometric features of those patterns crucially depend on the degree of dissociation and on the initial concentrations of the components. To a certain degree, the impact of the dissociation on the phase pattern formation effectively resembles a repulsion between the organic aggregates.
Finally, we emphasize that our computer simulations were performed within the framework of a very strongly simplified description of the reaction–diffusion processes. To further substantiate our predictions, a more sophisticated modeling is clearly needed. Especially, this concerns a more precise description of the nonlinear diffusion and phase separation of all the components by means of a more appropriate functional form of the Gibbs free energy and by including a concentration dependence of the diffusion coefficients. This is particularly important for larger polymer volume fractions and dissociation rates. Our current approach represents, thus, a first step aiming at demonstrating the potential of models combining phase separation and reaction–diffusion processes for elucidating the hypothesized basic mechanisms in diatoms morphogenesis. Moreover, additional experimental information concerning the in vivo conditions in the SDV will be needed to further fine tune our model towards more realistic conditions. Still, our approach can be of interest for studies aiming at inducing structure formation under in vitro conditions. The present simulations can also serve as a starting point for the modeling of silica structures, where the arrangement of organic aggregates provides a template for silica formation after addition of monosilicic acid to the system. We also hope to benefit from atomistic information obtained e.g. from Molecular Dyanmics simulations of the organic–inorganic interface properties in order to further improve the quality of the continuum model.
Availability of data and materials
Simulation data can be obtained upon request by the authors
Abbreviations
 SDV:

silica deposition vesicle
 LCPA:

longchain polyamine
References
 1.
Poulsen N, Kröger N. Silica morphogenesis by alternative processing of silaffins in the diatom Thalassiosira pseudonana. J Biol Chem. 2004;279(41):42993–9. https://doi.org/10.1074/jbc.M407734200.
 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. https://doi.org/10.1002/adfm.200500616.
 3.
Kröger N. Prescribing diatom morphology: toward genetic engineering of biological nanomaterials. Curr Opin Chem Biol. 2007;11(6):662–9.
 4.
Kröger N, Sandhage KH. From diatom biomolecules to bioinspired syntheses of silica and titaniabased materials. MRS Bull. 2010;35(2):122–6.
 5.
Tesson B, Hildebrand M. Characterization and localization of insoluble organic matrices associated with diatom cell walls: insight into their roles during cell wall formation. PLoS ONE. 2013;8(4):e61675. https://doi.org/10.1371/journal.pone.0061675.
 6.
Hildebrand M, Lerch SJL, Shrestha RP. Understanding diatom cell wall silicificationmoving forward. Front Mar Sci. 2018;5:1–19.
 7.
Pawolski D, Heintze C, Mey I, Steinem C, Kröger N. Reconstituting the formation of hierarchically porous silica patterns using diatom biomolecules. J Struct Biol. 2018;204(1):64–74. https://doi.org/10.1016/j.jsb.2018.07.005.
 8.
Bentley K, Cox EJ, Bentley PJ. Nature’s batik: a computer evolution model of diatom valve morphogenesis. J Nanosci Nanotechnol. 2005;5(1):25–34.
 9.
Gordon R, Drum RW. The chemical basis of diatom morphogenesis. Int Rev Cytol. 1994;150:243–372.
 10.
Lenoci L, Camp PJ. Diatom structures templated by phaseseparated fluids. Langmuir. 2008;24(1):217–23.
 11.
Parkinson J, Brechet Y, Gordon R. Centric diatom morphogenesis: a model based on a DLA algorithm investigating the potential role of microtubules. Biochim Biophys Acta Mol Cell Res. 1999;1452(1):89–102.
 12.
Willis L, Page KM, Broomhead DS, Cox EJ. Discrete freeboundary reaction–diffusion model of diatom pore occlusions. Plant Ecol Evol. 2010;143(3):297–306.
 13.
Willis L, Cox EJ, Duke T. A simple probabilistic model of submicroscopic diatom morphogenesis. J R Soc Interface. 2013;10(83):20130067.
 14.
Sumper M. A phase separation model for the nanopatterning of diatom biosilica. Science. 2002;295(5564):2430–3.
 15.
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.
 16.
Wieneke R, Bernecker A, Riedel R, Sumper M, Steinem C, Geyer A. Silica precipitation with synthetic silaffin peptides. Organ Biomol Chem. 2011;9(15):5482.
 17.
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.
 18.
Glotzer SC, Di Marzio EA, Muthukumar M. Reactioncontrolled morphology of phaseseparating mixtures. Phys Rev Lett. 1995;74(11):2034–7.
 19.
Kröger N, Lorenz S, Brunner E, Sumper M. Selfassembly of highly phosphorylated silaffins and their function in biosilica morphogenesis. Science. 2002;298(5593):584–6.
 20.
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.
 21.
Abacilar M, Daus F, Haas C, Brückner SI, Brunner E, Armin GA. Synthesis and NMR analysis of 13C and 15Nlabeled longchain polyamines (LCPAs). RSC Adv. 2016;6:93343–8.
 22.
Sumper M, Lorenz S, Brunner E. Biomimetic control of size in the polyaminedirected formation of silica nanospheres. Angew Chem Int Ed Engl. 2003;42(42):5192–5.
 23.
Lutz K, Gröger C, Manfred Sumper M, Brunner E. Biomimetic silica formation: analysis of the phosphateinduced selfassembly of polyamines. Phys Chem Chem Phys. 2005;7(14):2812–5.
 24.
Jantschke A, Koers E, Mance D, Weingarth M, Brunner E, Baldus M. Insight into the supramolecular architecture of intact diatom biosilica from DNPsupported solidstate NMR spectroscopy. Angew Chem Int Edn. 2015;54(50):15069–73.
 25.
Cahn JW. Free energy of a nonuniform system. II. Thermodynamic basis. J Chem Phys. 1959;30(5):1121–4.
 26.
Haasen P. Phase transformations in materials. Weinheim: VCH; 1990.
 27.
Ariyapadi MV, Nauman EB. Gradient energy parameters for polymer–polymersolvent systems and their application to spinodal decomposition in true ternary systems. J Polym Sci Part B Polym Phys. 1990;28(12):2395–409.
 28.
Zielinski JM, Duda JL. Predicting polymer/solvent diffusion coefficients using freevolume theory. AIChE J. 1992;38(3):405–15.
 29.
Masaro L, Zhu XX. Physical models of diffusion for polymer solutions, gels and solids. Progress Polym Sci. 1999;24(5):731–75.
 30.
Fares HM, Schlenoff JB. Diffusion of sites versus polymers in polyelectrolyte complexes and multilayers. J Am Chem Soc. 2017;139(41):14656–67.
 31.
Chen LQ, Shen J. Applications of semiimplicit Fourierspectral method to phase field equations. Comput Phys Commun. 1998;108(2–3):147–58.
 32.
https://www.comsol.de/. pp. 15069–15073.
 33.
Gräb O, Abacilar M, Daus F, Geyer A, Steinem C. 3Dmembrane stacks on supported membranes composed of diatom lipids induced by longchain polyamines. Langmuir. 2016;32(39):10144–52.
Acknowledgements
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 thank Dr. Andrea Benassi for his fruitful cooperation in an early stage of this study.
Funding
Deutsche Forschungsgemeinschaft: Nanopatterned Organic Matrices in Biological Silica Mineralization Research Unit FOR 2038). The DFG only provided funding for the positions of H. E.. 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. Open access funding provided by Projekt DEAL.
Author information
Affiliations
Contributions
MB, HY, and DW developed the computational model and implemented the simulations. RG, AD, HE and GC wrote and corrected the manuscript. All authors have read and approved the manuscript.
Corresponding author
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.
2D concentration simulation with permanent coarsening DOI: 10.6084/m9.figshare.12021081 (still private). 2D concentration profile simulation of component C for a parameter set where permanent coarsening of phase regions occurred. (see Fig. 4). 2D concentration simulation with stable pattern DOI: 10.6084/m9.figshare.12021033 (still private). Simulation of a 2D concentration profile of component C. High initial concentrations of A and B. (see Fig. 6b). Rough estimation of forward and backward reaction rates used in the manuscript and additional simulations with different parameter sets
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 http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Bobeth, M., Dianat, A., Gutierrez, R. et al. Continuum modelling of structure formation of biosilica patterns in diatoms. BMC Mat 2, 12 (2020). https://doi.org/10.1186/s42833020000187
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s42833020000187
Keywords
 Diatoms
 Biosilica formation
 Phase separation
 Cahn–Hilliard equation
 Reaction–diffusion equations