Continuum modelling of structure formation of biosilica patterns in diatoms

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. Phase-separation of species-specific 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 self-assembly 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 three-dimensional 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 diffusion-limited 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 diffusion-reaction mechanisms, the formation of various nano-sized 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 self-assembly 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 two-dimensional 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 pre-existing silica costae within the base layer. Depending on the choice of this regular pre-patterning field, different highly symmetric valve structures could be achieved.

Extensive experimental investigations have revealed that species-specific 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, filament-like 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 network-like silica morphology was formed. The long-term 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 far-reaching 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 self-assembly of long-chain polyamines (LCPAs). Besides the experimental studies mentioned above, the importance of phosphate ions is further highlighted by experiments on silaffins, which showed that silaffin 1-A (in contrast to natSil-1A) was not able to promote silica precipitation in the absence of phosphates [19]. This was traced back to the absence in silaffin 1-A 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:

$$\begin{aligned} \hbox {A} + \hbox {B} \rightleftharpoons \hbox {C}, \end{aligned}$$

(1)

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 :

$$\begin{aligned} \hbox {H}_{2}\hbox {PO}^{-}_{4}+\hbox {C}_{10}\hbox {H}_{27}\hbox {N}_{6}^{\rm{+n}}\leftrightarrow [\hbox {H}_{2}\hbox {PO}^{-}_{4}+ \hbox {C}_{10}\hbox {H}_{27}\hbox {N}_{6}]^{+(\rm{n-1})}. \end{aligned}$$

To describe the phase separation process, we use the well-known Cahn-Hilliard model [10, 25, 26]. We start from the Gibbs energy functional:

$$\begin{aligned} G[\phi ]=\int {\rm{d}}V \left[ g(\phi )+\frac{1}{2}\kappa (\nabla \phi )^2\right] , \end{aligned}$$

(2)

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

$$\begin{aligned} \dfrac{\partial \phi }{\partial {t}}&= \nabla \cdot M\nabla \mu , \end{aligned}$$

(3)

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]:

$$\begin{aligned} \dfrac{\partial \phi }{\partial {t}}&= \nabla \cdot M\nabla \left( g'(\phi )-\kappa \nabla ^2\phi \right) . \end{aligned}$$

(4)

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 Flory-Huggins model:

$$\begin{aligned} g(\phi )&=g_0\left[ \frac{1}{N}\phi \ln \phi +(1-\phi )\ln (1-\phi )\right. \nonumber \\&\qquad \left. \quad +\chi \phi (1-\phi )\right] , \end{aligned}$$

(5)

where $g_0$ is a prefactor, $N$ is the degree of polymerization, and $\chi$ is the Flory-Huggins 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:

$$\begin{aligned} \dfrac{\partial \phi }{\partial {\hat{t}}}&= \hat{\nabla }^2\left( \hat{g}'(\phi )-\hat{\nabla }^2\phi \right) , \end{aligned}$$

(6)

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]:

$$\begin{aligned} \dfrac{\partial c_{\rm{A}}}{\partial {t}}&= \nabla \cdot D_{\rm{A}}\nabla c_{\rm{A}}-\alpha c_{\rm{A}} c_{\rm{B}} + \beta c_{\rm{C}}, \end{aligned}$$

(7)

$$\begin{aligned} \dfrac{\partial c_{\rm{B}}}{\partial {t}}&= \nabla \cdot D_{\rm{B}}\nabla c_{\rm{B}}-\alpha c_{\rm{A}} c_{\rm{B}} + \beta c_{\rm{C}}, \end{aligned}$$

(8)

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 polyamine-phosphate ion-water system, we will assume in the following that these kinetic coefficients are constants in a first analysis. A rough order-of-magnitude 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:

$$\begin{aligned} \dfrac{\partial \hat{c}_{\rm{A}}}{\partial {\hat{t}}}&= \hat{\nabla }\cdot \hat{D}_{\rm{A}}\hat{\nabla }^2 \hat{c}_{\rm{A}}-\hat{\alpha } \hat{c}_{\rm{A}} \hat{c}_{\rm{B}} + \hat{\beta } \phi \end{aligned}$$

(9)

$$\begin{aligned} \dfrac{\partial \hat{c}_{\rm{B}}}{\partial {\hat{t}}}&= \hat{\nabla }\cdot \hat{D}_{\rm{B}}\hat{\nabla }^2 \hat{c}_{\rm{B}}-\hat{\alpha } \hat{c}_{\rm{A}} \hat{c}_{\rm{B}} + \hat{\beta } \phi \end{aligned}$$

(10)

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:

$$\begin{aligned} \dfrac{\partial \phi }{\partial {\hat{t}}} = \hat{\nabla }^2\left( \hat{g}'(\phi )-\hat{\nabla }^2\phi \right) +\hat{\alpha } \hat{c}_{\rm{A}} \hat{c}_{\rm{B}} - \hat{\beta } \phi . \end{aligned}$$

(11)

For the sake of brevity, we will omit the roof on the scaled variables in the following.

Fig. 1

Dispersion relation $\omega (q)$ a for different values of $\beta = 0, 0.05, 0.1, 0.2$ and $c_{\rm{B,0}} = 0.25$ revealing the strong impact of the dissociation constant on the spinodal decomposition; b for different initial concentration of component B (phosphate ions) $c_{\rm{B,0}} = 0.1, 0.15, 0.2, 0.4$ and $\beta = 0.1$. Shared parameters: $D_{\rm{A}}={100}$, $D_{\rm{B}}={1000}$, $\alpha =10$ , $c_{\rm{A,0}} = 0.25$

It should be noted that the chosen expression for the Gibbs energy in Eq. 5 is valid for a polymer-solvent 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:

$$\begin{aligned} \phi _e&= \frac{1}{2}\left[ c_{\rm{A,0}} + c_{\rm{B,0}} + \frac{\beta }{\alpha } \right. \nonumber \\&\quad \left. - \sqrt{\left( c_{\rm{A,0}} +c_{\rm{B,0}} + \frac{\beta }{\alpha } \right) ^2 - 4c_{\rm{A,0}}c_{\rm{B,0}}}\, \right] . \end{aligned}$$

(12)

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 semi-implicit Fourier-spectral 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

Fig. 2

1D concentration profiles: Two snapshots of concentration profiles of components A, B, and C at two times. Parameters: $D_{\rm{A}}=2$, $D_{\rm{B}}=10$, $\alpha =10$, $c_{\rm{A,0}} =c_{\rm{B,0}} = 0.15$, $\beta = 0.02$

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:

$$\begin{aligned} \frac{\partial \delta c_{\rm{A}}}{\partial t}&= D_{\rm{A}}\nabla ^2 \delta c_{\rm{A}}\nonumber \\&\quad -\alpha \left( \delta c_{\rm{A}} c_{\rm{Be}} + c_{\rm{Ae}} \delta c_{\rm{B}}\right) +\beta \delta \phi \end{aligned}$$

(13)

$$\begin{aligned} \frac{\partial \delta c_{\rm{B}}}{\partial t}&= D_{\rm{B}}\nabla ^2 \delta c_{\rm{B}}\nonumber \\&\quad -\alpha \left( \delta c_{\rm{A}} c_{\rm{Be}} + c_{\rm{Ae}} \delta c_{\rm{B}}\right) +\beta \delta \phi \end{aligned}$$

(14)

$$\begin{aligned} \dfrac{\partial \delta \phi }{\partial {t}}&= \nabla ^2\left( g''(\phi _{\rm{e}})\delta \phi -\nabla ^2\delta \phi \right) \nonumber \\&\quad +\alpha \left( \delta c_{\rm{A}} c_{\rm{Be}} + c_{\rm{Ae}} \delta c_{\rm{B}}\right) - \beta \delta \phi . \end{aligned}$$

(15)

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:

$$ \begin{array}{ll} \left( {\begin{array}{*{20}c} {m_{{\text{A}}} } & {\alpha c_{{{\text{Ae}}}} } & { - \beta } \\ {\alpha c_{{{\text{Be}}}} } & {m_{{\text{B}}} } & { - \beta } \\ { - \alpha c_{{{\text{Be}}}} } & { - \alpha c_{{{\text{Ae}}}} } & {m_{\phi } } \\ \end{array} } \right)\left( {\begin{array}{*{20}c} {a_{{\text{A}}} } \\ {a_{{\text{B}}} } \\ {a_{\phi } } \\ \end{array} } \right) = 0 \hfill \\ m_{{\text{A}}} (q) = \omega + D_{{\text{A}}} q^{2} + \alpha c_{{{\text{Be}}}} \hfill \\ m_{{\text{B}}} (q) = \omega + D_{{\text{B}}} q^{2} + \alpha c_{{{\text{Ae}}}} \hfill \\ m_{\phi } (q) = \omega + q^{2} \left( {g^{{\prime \prime }} (\phi _{{\text{e}}} ) + q^{2} } \right) + \beta . \hfill \\ \end{array} $$

(16)

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 q-values 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 long-wavelength 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

Fig. 3

1D concentration profiles of component C for different dissociation constants $\beta = 0.15$ and $\beta = 0$. Without dissociation considerable coarsening of phase regions occurred (initial concentrations $c_{\rm{A,0}} =c_{\rm{B,0}} = 0.35$)

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.

Fig. 4

2D concentration profile snapshots of component C for a parameter set where permanent coarsening of phase regions occurred $c_{\rm{A,0}} = c_{\rm{B,0}} =0.2$, $\beta = 10^{-4}$. For animation see Additional file 1

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 one-dimensional 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, two-dimensional (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.

Fig. 5

Stable 2D concentration profile: Stationary regular pattern of concentration C ($c_{\rm{A,0}} = c_{\rm{B,0}} =0.2$, $\beta = 0.1$). Note the stability of the pattern up to a very large computation time

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 organic-rich phase of component C with circular pores of organic-poor phase develops. The size of the pores increases with decreasing dissociation constant. Interestingly, at certain values of the dissociation constant, also lamellar (labyrinth-like) phase patterns were obtained. It is remarkable that those stripes are stationary despite their irregular structure.

Fig. 6

Stationary 2D concentration profiles of component C at higher initial concentrations of A and B for different values of $\beta$. Parameters: a $c_{\rm{A,0}} = c_{\rm{B,0}} =0.4$, $t=274259$, b $c_{\rm{A,0}} = c_{\rm{B,0}} =0.4$, $t=535837$, c) $c_{\rm{A,0}}=c_{\rm{B,0}}=0.5$, $t=43385$. For animation of case b) see Additional file 1

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], pre-existing large pores represent confinements for the formation of smaller pores. Pattern formation in pre-existing 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, organic-rich 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.

Fig. 7

Circular stationary patterns: Concentration profiles of component C within a circular boundary for different dissociation constants (columns, $\beta = 0.15, 0.1, 0.01, 0.001$) and initial concentrations of components A and B ( $c_{\rm{A,0}} = c_{\rm{B,0}} = 0.2, 0.25, 0.4, 0.5$). The radius of the simulation cell was $r = 40$. Bright regions correspond to the organic-rich phase

Fig. 8

Attractive interaction: Snapshots of aggregation of component C (time $t = 140, 200, 230, 500$) influenced by an attractive interaction of component A with the left boundary, which leads to the formation of an oriented hexagonal pattern. Parameters: $h_0 = -0.01$, $c_{\rm{A,0}} = c_{\rm{B,0}} = 0.2$, $\beta = 0.1$ , $\xi = 1$

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 long-chain 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

Fig. 9

3d patterns: Iso-concentration surfaces and corresponding surface plots of the concentration of component C in a square based box, showing the formation of stationary cylindrical in a as well as stationary half-sphere aggregates at the top and bottom face of the simulation cell in b. The thickness of the simulation cell in b of $\hbox {d} = 7$ is slightly thinner than the thickness of $\hbox {d} = 8$ in a. Parameters: $c_{\rm{A,0}} = c_{\rm{B,0}} = 0.25$, $\beta =0.15$ and $t=5 \cdot 10^{5}$

Microscopic images of the three-dimensional 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 z-direction perpendicular to the flat SDV. Periodic boundary conditions were supposed only in x- and y-directions, whereas at the boundaries in z-direction 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 iso-concentration 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 half-sphere 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 pre-patterning 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.