research papers
Crystallographic phase retrieval method for
bicontinuous phases: indicator-based methodaDepartment of Physics, Faculty of Science, Shizuoka University, Shizuoka, 422-8529, Japan, and bNanomaterials Research Division, Research Institute of Electronics, Shizuoka University, Shizuoka, 422-8529, Japan
*Correspondence e-mail: oka.toshihiko@shizuoka.ac.jp
An indicator-based crystallographic phase retrieval method has been developed for diffraction data of bicontinuous cubic phases of lyotropic liquid crystals. Such liquid crystals have large structural disorder; the number of independent Bragg reflections that can be observed is limited. This paper proposes two indicators to identify plausible combination(s) of crystallographic phases, i.e. electron-density distribution. The indicators are based on the characteristics of the liquid crystals: molecules diffuse mainly in the direction parallel to polar–nonpolar interfaces and the electron density in the direction parallel to the interfaces is averaged. One indicator is the difference between the maximum and minimum electron density, and the other is calculated from the Hessian matrix of the electron density. Using test data, the electron densities were calculated for all possible phase combinations and indicators were obtained. The results indicated that the electron densities with the minimum indicators were close to the true electron density. Therefore, this method is effective for phase retrieval. The accuracy of the phase retrieval decreased when the of the region including the triply periodic minimal surface increased.
Keywords: lyotropic liquid crystals; triply periodic minimal surfaces; crystallographic phase retrieval.
1. Introduction
Restoring the phases of the structure factors is a crucial problem in ; Giacovazzo, 2001). Closely related is the charge-flipping method, an efficient iterative phase retrieval method with similar requirements that has become popular in recent years (Oszlányi & Sütő, 2004, 2008). It may be possible to obtain the structure using data from samples that satisfy the structural information; however, this method is limited to solid crystals in which the atomic positions are almost fixed.
determination from X-ray diffraction (XRD) data. solve the by combining the intensity information with constraints based on the expected structural information (atomicity, positivity and zero density regions) (Sayre, 1952Liquid crystals have properties that are intermediate between solids and liquids. In liquid crystals, the molecules move like a liquid in a certain direction, but are oriented like a solid crystal in other directions. For the latter reason, many et al., 1996), and diffraction data can be obtained as if they were 3D crystals [Fig. 1(a)]. However, the spatial resolution of the observed XRD data is low due to structural disorder, and the number of independent reflections is small. The previously mentioned constraints, except the positivity, are not satisfied. There is thus a problem in restoring the electron density from XRD data, although the centrosymmetry of the LLC bicontinuous phase restricts the phase to 0 or π. Luzzatti and colleagues first clarified the electron-density distribution of an LLC bicontinuous phase (Luzzati et al., 1988; Mariani et al., 1988). They proposed to use the average of the fourth moment of the electron density as an indicator for phase retrieval. They insisted that the combination of the phases at the smallest value was the most reliable. However, as they admitted, the electron density with the minimum indicator was not always the proper phase combination. Other methods have also been proposed for LLCs, such as looking for the position of the methyl group in the hydrocarbon chain of a lipid bilayer (Harper et al., 2000).
phases have periodic structures and are subject to XRD measurements. Bicontinuous cubic phases of lyotropic liquid crystals (LLCs) have structures similar to triply periodic minimal surfaces (TPMSs) (HydeWe have performed X-ray crystallographic studies of LLC bicontinuous cubic phases using single-crystal regions of samples (Oka, 2017; Oka et al., 2018, 2020). In these studies, the structures of the bicontinuous cubic phases were determined using models based on available information. However, it would be better if the structure of the bicontinuous cubic phase could be revealed without a structural model. Here, a simple indicator-based phase retrieval method is reported. Two valuable indicators have been identified based on the universal characteristics of the LLC bicontinuous structures. The phase combination in which these indicators are minimized is very close to the true phase combination, and makes it possible to retrieve the crystallographic phase of the TPMS-like structure. Although this method was designed to clarify the structure of the LLC bicontinuous cubic phases, it would be applicable to thermotropic liquid crystals, polymers, and other materials with TPMS-like structures.
2. Indicators for bicontinuous structure (TPMS-like structure)
2.1. Bicontinuous structure (TPMS-like structure)
An LLC consists of two or more components, ). In type II bicontinuous cubic phases of LLCs, three different TPMS-like structures, P surface, D surface and gyroid (G) surface, have been observed, while only the G surface has been observed in type I (Hyde et al., 1996). The space groups of the three TPMS-like structures of P, D and G are , and , respectively. In type II, the nonpolar region contains the TPMS as the center surface, and the polar region is located on the two networks separated by the TPMS [Fig. 1(a)]. In the type I bicontinuous cubic phase, the positions of the polar and nonpolar regions are swapped with those of type II (Israelachvili, 2011). In both cases, the polar–nonpolar interface is on the molecule, and the molecule can move like a liquid molecule in the direction parallel to the polar–nonpolar interface. Therefore, the electron density in the direction parallel to the interface is averaged, but not that in the perpendicular direction. This averaging leads to the separation of the polar and nonpolar regions. In the case of an molecule and water, if the electron densities of the polar parts of the molecules and waters are almost the same, then the electron density can be described in two levels, polar and nonpolar, without considering fluctuations. The region in the TPMS is a minimal surface like structure, and the regions on the networks are cylinder-like structures with branches (Fig. 2). If the electron densities of the polar parts of water and amphiphiles are different, then the electron density can be described in three continuous regions, located on or along the TPMS or the networks.
molecules and water (or oil). There are two types of LLC, oil-in-water and water-in-oil types, commonly referred to as type I and type II, respectively (Israelachvili, 20112.2. Indicator Iρ
The bicontinuous structure consists of two (or three) flat electron-density regions if fluctuations are not considered [Fig. 1(b)]. If electron density is restored using structure factors with incorrect phases, then the density in the flat region will undulate and the difference between the minimum and maximum electron densities will become more significant compared with the true density. Therefore, the difference between the minimum electron density in a , and the maximum electron density, , would serve as an indicator:
This indicator corresponds to the range in statistics. A smaller indicator is a better candidate for the proper phase solution. The actual electron density is close to the step-like electron density with the convolution of the fluctuation function [Fig. 1(c)]. In this case, the electron density is almost constant at a maximum or minimum on the TPMS and the networks [Fig. 1(a)]. Electron densities calculated with incorrect phases will undulate on the TPMS or networks, which increases the indicator .
2.3. Indicator IK
A minimal surface like structure located on the TPMS and cylinder-like structures with junctions located on the networks are close to parts of the electron density of the LLC bicontinuous cubic phases [Figs. 2(b) and 2(c)]. The molecules diffuse along the structural motifs; therefore, the isoelectron density surfaces are not closed but infinitely connected. However, only the sphere-like structure in Fig. 2 has closed isoelectron density surfaces, i.e. a region that is strictly convex upward. An isolated density maximum/minimum is difficult to imagine in the LLC bicontinuous cubic phase because it seemingly contradicts molecular diffusion over a long distance. Therefore, few density regions that are convex upward or downward are expected. Convex regions can be determined by the eigenvalues of the Hessian matrix of the electron density, :
where subscripts indicate partial derivatives. If all eigenvalues of the Hessian matrix are positive, then the region is strictly convex downward, and if all are negative, then the region is strictly convex upward (Rockafellar & Wets, 2010). Let C be the electron-density regions that are convex upward or downward and let K(r) be the determinant of H(r): K(r) = det[H(r)]. Therefore, K(r) is the product of three eigenvalues. The following indicator was considered:
which is the integral of the absolute value of in the region C of the electron density in which the three eigenvalues are of the same sign in a where r is the position in a Thus, IK is an indicator that expresses the total (integrated) convexity. With regard to the geometrical meaning, the eigenvalues of the Hessian matrix are the principal curvatures, and K(r) is the Gaussian curvature if the 3D electron density is regarded as a hypersurface in a 4D hyperspace (Monga & Benayoun, 1995; Goldman, 2005).
3. Phase retrieval method based on the Iρ and IK indicators
The
can be divided into the amplitude and the phase part , where is the phase. In a diffraction experiment, intensities are measured as the squares of the amplitudes and the phases are not observable. The electron density can be obtained by Fourier transformation of the structure factors as follows:where r is the position in the h is the reciprocal-lattice vector and V is the volume of the Hobs is the set of h for which the could be observed. F(000) was not used in this paper; therefore, the integrated value of the electron density in the is zero. High-electron-density regions are thus positive, whereas low-electron-density regions are expressed as negative in .
In the phase retrieval, the electron densities were calculated for all possible phase combinations by placing them into the phase part of the π (if the unit-cell origin is placed on a symmetry center), i.e. the structure factors are real (positive or negative, respectively). When n is the number of the independent structure-factor amplitude, the corresponding number of phase combinations is 2n. The electron density was calculated in 32 × 32 × 32 voxels as a was obtained from the difference between the maximum and minimum electron densities. To obtain IK, the Hessian matrices were obtained from the electron densities for all voxel points in the which yielded three eigenvalues at all points. The region where the signs of the three eigenvalues were the same was designated as C, and K(r) was obtained at each point within C. IK was then calculated by integration as described in Section 2.3. The fourth moment of the electron density was calculated at the same time (Luzzati et al., 1988). F(000) was not used in this paper; therefore, corresponds to the average of the fourth power of the electron density, and was calculated as the average of the fourth power of the electron density in the unit cell.
The structures of the bicontinuous cubic phases are centrosymmetric, so that the phase values are 0 orAccording to the Babinet principle, when the electron density is inverted, the amplitudes of the structure factors are the same, but the phases are shifted by π. Therefore, the structures with inverted electron densities were considered equivalent, and the number of phase combinations to be calculated was reduced. In the , if the origin of a structure is (0,0,0), then moving the origin to (1/2,1/2,1/2) will give equivalent amplitudes of the structure factors, although some phases are shifted by π (Giacovazzo, 2001). Therefore, it was considered that the structure is equivalent, even when the translation operation of (1/2,1/2,1/2) is performed, and the number of phase combinations to be calculated is reduced. The electron densities and indicators were calculated using these structure factors. Phase retrieval was conducted using an in-house-made Python3 script.
4. Test data and quality of phase retrieval
4.1. Test data derived from constructed models
Electron-density models were constructed to test the indicators for the phase retrieval. Two polar–nonpolar interface models for the LLC bicontinuous cubic phase were used: a parallel surface (PS) (Hyde et al., 1996) and a constant mean curvature surface (CMCS) (Anderson et al., 1990; Große-Brauckmann, 1997). In the PS model, the polar–nonpolar interface is parallel to the TPMS, and in the CMCS model, the interface is a CMCS. The PS and CMCS models were generated as model electron densities from the three TPMSs, i.e. the G, P and D surfaces. The PS and CMCS models based on the G surface are referred to as the G-PS and G-CMCS models, respectively.
The TPMSs and CMCSs were created using Surface Evolver (Brakke, 1992; Shearman et al., 2007). The PSs were generated as surfaces at fixed distances from the TPMSs [Fig. 1(b)]. Stepwise electron densities were created for both models. The electron density was calculated in 128 × 128 × 128 voxels as a The density from the interface to the TPMS side was set to 1.0, and the rest of the region was set to 0.0. The volume fractions from the interface to the TPMS side () were 0.2, 0.4, 0.6, 0.7 and 0.8. PS models with different were generated by varying the constant distance from the TPMS to the polar–nonpolar interface (Qiu & Caffrey, 1998), while the CMCS model requires only a to generate the CMCS. The CMCS with could not be generated for the P surface; therefore, only P-CMCS models with 0.2 and 0.4 were used. Using values from previous papers as a reference (Oka, 2017; Oka et al., 2018, 2020), the lattice constant was set to 1 and the width of the Gaussian function was set to 0.05 for the G surface based models and 0.07 for the other models.
The structure factors were obtained by Fourier transformation of the model densities :
Integration is performed in a Mathematica 12.1 (Wolfram Research, Inc., USA).
The numbers of independent structure-factor amplitudes used were 22, 19 and 13 for the G, D and P surface based models, respectively. The model density construction, except for the formation of the TPMS and CMCS, and other calculations were conducted using4.2. Rp value to evaluate the quality of phase retrieval
The Rp value is defined as a quantity that evaluates how well the phase combination used in the calculation agrees with the true phase combination:
Rp close to 0 indicates a good agreement with the true phase combination. From the Babinet principle, the inverted electron density was also treated in the same way, and the π-shifted phase combination from was also treated as the true phase. Of the Rp calculated from the original and π-shifted phase combinations, that with the smaller value was adopted.
5. Results
5.1. Phase retrieval for constructed models
Phase retrieval was performed for model densities with known true phase combinations. The electron density was calculated for all phase combinations, and the corresponding IK and indicators were calculated.
Fig. 3(a) shows the distribution of IK and for all phase combinations at of the G-PS model. The points are distributed from the lower left to upper right. Therefore, there is an approximate positive correlation between IK and . The Rp values of the points tend to be smaller in the lower left region; points with Rp of approximately 0.1 or less are concentrated in the lower left region. IK and for the true phase combination are in the lower left corner, which indicates that both indicators of the true phase are almost the smallest, as expected; the same is true when is 0.2 and 0.6 in the other models (Fig. S1 in the supporting information). However, when is 0.8, IK for the true phase combination is roughly minimized, whereas is clearly not minimized [Fig. 3(b)]. The tendency that the indicator for the true phase combination is not minimized when is large is also true for the G-CMCS model and the other surface-based models (Fig. S1). For example, in the G-CMCS model, both IK and are not minimized when is 0.7 or 0.8.
The Rp values of the phase combinations with the minimum indicators in the G-PS and G-CMCS models are shown in Fig. 4. When the is smaller than 0.6, Rp of the phase combinations at the minimum indicators is approximately 0.1 or less. Therefore, in this range of volume fractions, both IK and indicators are useful for phase retrieval. When Vfrac is 0.7 or 0.8, some of the phase combinations at the minimum indicator are Rp > 0.2, which indicates that the phase retrieval by the indicator is not realized. In particular, the phase combination with the minimum differs significantly from the true combination for with all the models. This indicates that when Vfrac is large, the constraint that the difference between the minimum and maximum electron density is minimal does not hold. On the other hand, the Rp value for minimum IK is relatively small, even when the is 0.7 or 0.8. Including other models, Rp > 0.2 with minimum IK is limited to three conditions: the G-CMCS model with and the P-PS model with = 0.7 and 0.8. Therefore, IK is useful for phase retrieval to some extent, even if . When Vfrac is large in the G-CMCS model, the interface structure deviates from the simple model described in Section 2 due to the junctions of the networks (Große-Brauckmann, 1997). In the P-PS model, the junction has six branches, so if the is large, the junction is far outside the cylinder structure. The combination of these factors and the smearing caused by fluctuations leads to a failure in the phase retrieval of IK because the electron density is out of the constraint based on the eigenvalues.
The previously used indicator, the fourth moment of the electron density (Luzzati et al., 1988), was compared with . In all models, when Vfrac = 0.2, the Rp value of the minimum phase combination is larger than that of (Figs. 3 and S2). When Vfrac = 0.7 for D-PS and D-CMCS and Vfrac = 0.6 for P-PS, the Rp value for the minimum is larger than that based on . Although there are conditions in which the Rp values are similar for both indicators, gives better results than under several conditions. Both and are indicators based on electron density: is calculated from the entire electron density, while is calculated from the maximum and minimum electron densities. Therefore, uses information from the intermediate-electron-density region, which is not used in . This may cause the phase retrieval results of to be better under several conditions.
5.2. Phase retrieval for experimental data
The experimental data of LLC bicontinuous cubic phases acquired in previous X-ray single-crystal structure analyses (Oka, 2017; Oka et al., 2018, 2020) were used as test data. In the calculation of the Rp value, the phase combinations previously derived from the model analysis were taken as the true combinations.
Monoolein, a type of lipid, and water form three type II bicontinuous cubic phases with TPMSs of G, D and P (Oka, 2017). The number of independent structure factors is from 8 to 14. can be estimated to be around 0.43 to 0.54 (Table S2). The scatter plot shows that IK and are roughly positively correlated, and that the Rp value is smaller toward the lower left [Fig. S3(a)]. This is the same as when in the model densities. The phase combinations obtained from the structural models in the previous paper are also located at the lower left of the scatter plot. In the phase combinations with the minimum indicators, the calculated Rp values were almost zero (Table 1). This shows a good agreement between the structure shown by the indicator-based phase retrieval method and the structural model in the previous paper (Oka, 2017). For the two type II bicontinuous cubic phases of phytantriol with TPMSs of G and D, the number of independent reflections is 21, and the values are 0.57 and 0.66 (Oka et al., 2018). These data are also similar to those of monoolein, and the structure shown by the phase retrieval method and the structure model in the previous paper (Oka et al., 2018) are in good agreement (Table 1).
Hexaethylene glycol monododecyl ether (C12EO6) and water form a type I bicontinuous cubic phase of G TPMS, with (Oka et al., 2020). In the scatter plot of IK and , the distribution of points is split in two on the lower left side [Fig. S3(c)]. IK is a minimum near the point of the phase combination that was considered true in the previous paper (Oka et al., 2020) (), whereas is not a minimum near the point. On the other hand, close to the point where is minimum (), Rp is approximately 0.8 to 0.9. The corresponding electron densities are shown in Fig. 5. The electron density of the minimum IK is in very good agreement with the electron density in the previous paper (Oka et al., 2020). A comparison of the electron density of the minimum IK [Fig. 5(b)] with Figs. 1(b) and 1(c) indicates that the width of the high-electron-density region on the TPMS side is wider, which is consistent with the large on the TPMS side of C12EO6. With the minimum , the high-electron-density region on the TPMS side is narrower and the electron-density regions on the network sides are larger compared with the minimum IK. If the is known in advance, then the electron density indicates that the phase combination with the minimum is not close to the true phase combination.
The fourth moment of the electron density was also calculated in the experimental data, and the results are summarized in Table 1. The Rp values for the phase combination with the minimum indicator are comparable with those of for four of the six experimental data sets. On the other hand, the Rp value of the minimum indicator is approximately 0.1 for monoolein P and 0.2 for D, which are larger than that of . For the difference of Rp ∼ 0.1, the electron densities show similar structures, but there are non-negligible differences in the details (Fig. S6). Considering the results with the model structures, IK is the best indicator with the widest range of applicable conditions, followed by and then . Regardless of which indicator is used, care should be taken because there are volume fractions for which phase retrieval is difficult.
Unobservable structure factors cause truncation errors in the Fourier synthesis of electron density, . The truncation errors may make the indicators inaccurate. However, given the observed amplitudes, the amplitudes of unobservable structure factors are probably less than 1% of the maximum amplitude. The amplitude decreases rapidly at high scattering angles due to the structural disorder of the i.e. indicators, seems insignificant. On the other hand, it may affect the phase retrieval of structure factors with weak amplitudes and cause the Rp value to deviate slightly from 0.
(this is why the number of observable amplitudes is small). For this reason, the effect of the truncation error on the resultant structure model,6. Conclusion
The indicators and IK were determined to be useful for the crystallographic phase retrieval in LLC bicontinuous phases. The phase combination with the minimum indicator was very close to the true phase combination when the was less than 0.6. Even for volume fractions of approximately 0.7 to 0.8, IK was effective to determine the proper phase combination as an indicator, although there were exceptions. In the scatter plots of IK and , when is small, the points extending to the lower left do not spread much [Fig. 3(a)]. On the other hand, when is large, the points in the lower left corner branch and spread out [Fig. 3(b)], which seems to be related to the certainty of the indicator.
The model used for the phase retrieval was a simple structure with different electron densities on the TPMS side and the network side. On the other hand, the sample of experimental data used for the phase retrieval has a more complicated structure. In the case of a type I bicontinuous cubic phase such as C12EO6, the polar region is located on the TPMS side and the nonpolar region is located on the network side. The hydrophobic chain of the molecule is located in the nonpolar region, water is located in the TPMS in the polar region, and the hydrophilic part of the molecule is located in the polar–nonpolar interface side of the polar region. Therefore, the electron density differs among the three regions, with the hydrophilic part of the molecule having the highest electron density [Fig. 5(b)]. Phase retrieval was possible even in such a system.
It takes a long time to calculate the indicators of electron densities for all phase combinations; especially for IK, it was necessary to calculate the eigenvalues of the Hessian matrix at each voxel, which takes a long time. On the other hand, the structures of the LLC bicontinuous cubic phase have highly symmetric space groups of , and . Therefore, the computational time can be reduced by calculating only in the asymmetric units. In addition, the Babinet principle does not allow positive and negative electron densities to be distinguished in this paper; therefore, even when they are reversed, the number of combinations can be reduced. For an LLC bicontinuous cubic phase, the largest number of independent structure factors so far is 21 for phytantriol and C12EO6–8 (Oka et al., 2018, 2020). When the number of voxels was 323, the computation times of the Python script were 16 to 20 min using an AMD 3950X CPU with 16 cores and 32 threads. Calculation of 230 phase combinations under the same conditions would thus take 10 to 14 days. It may be possible to increase the calculation speed by revising the algorithm. If the IK calculation is eliminated, then the calculation time will be much shorter.
Phase retrieval was performed using models and experimental data of LLC bicontinuous cubic phases. TPMS-like structures have been observed in thermotropic liquid crystals, polymers and other systems; it should also be possible to apply this phase retrieval method to these systems.
Supporting information
Supporting information. DOI: https://doi.org/10.1107/S2053273322006970/ik5003sup1.pdf
Funding information
The following funding is acknowledged: JSPS KAKENHI (grant Nos. JP20H04634, JP18K03557).
References
Anderson, D. M., Davis, H. T., Scriven, L. E. & Nitsche, J. C. C. (1990). Advances in Chemical Physics, edited by I. Prigogine & S. A. Rice, pp. 337–396. New York: John Wiley & Sons, Inc. Google Scholar
Brakke, K. A. (1992). Exp. Math. 1, 141–165. CrossRef Google Scholar
Giacovazzo, C. (2001). International Tables for Crystallography, Vol. B, Reciprocal Space, edited by U. Shmueli, pp. 210–234. Dordrecht: Springer Netherlands. Google Scholar
Goldman, R. (2005). Comput. Aided Geom. Des. 22, 632–658. Web of Science CrossRef Google Scholar
Große-Brauckmann, K. (1997). Exp. Math. 6, 33–50. Google Scholar
Harper, P. E., Gruner, S. M., Lewis, R. N. A. H. & McElhaney, R. N. (2000). Eur. Phys. J. E, 2, 229–245. Web of Science CrossRef CAS Google Scholar
Hyde, S., Blum, Z., Landh, T., Lidin, S., Ninham, B. W., Andersson, S. & Larsson, K. (1996). The Language of Shape: the Role of Curvature in Condensed Matter: Physics, Chemistry and Biology. Amsterdam: Elsevier. Google Scholar
Israelachvili, J. N. (2011). Intermolecular and Surface Forces, revised 3rd ed. Boston: Academic Press. Google Scholar
Luzzati, V., Mariani, P. & Delacroix, H. (1988). Makromol. Chem. Macromol. Symp. 15, 1–17. CrossRef Web of Science Google Scholar
Mariani, P., Luzzati, V. & Delacroix, H. (1988). J. Mol. Biol. 204, 165–189. CrossRef CAS PubMed Web of Science Google Scholar
Monga, O. & Benayoun, S. (1995). Comput. Vis. Image Underst. 61, 171–189. CrossRef Web of Science Google Scholar
Oka, T. (2017). J. Phys. Chem. B, 121, 11399–11409. Web of Science CrossRef CAS PubMed Google Scholar
Oka, T., Ohta, N. & Hyde, S. (2018). Langmuir, 34, 15462–15469. Web of Science CrossRef CAS PubMed Google Scholar
Oka, T., Ohta, N. & Hyde, S. T. (2020). Langmuir, 36, 8687–8694. Web of Science CrossRef CAS PubMed Google Scholar
Oszlányi, G. & Sütő, A. (2004). Acta Cryst. A60, 134–141. Web of Science CrossRef IUCr Journals Google Scholar
Oszlányi, G. & Sütő, A. (2008). Acta Cryst. A64, 123–134. Web of Science CrossRef IUCr Journals Google Scholar
Qiu, H. & Caffrey, M. (1998). J. Phys. Chem. B, 102, 4819–4829. Web of Science CrossRef CAS Google Scholar
Rockafellar, R. T. & Wets, R. J.-B. (2010). Variational Analysis. Heidelberg, Berlin: Springer. Google Scholar
Sayre, D. (1952). Acta Cryst. 5, 60–65. CrossRef CAS IUCr Journals Web of Science Google Scholar
Shearman, G. C., Khoo, B. J., Motherwell, M.-L., Brakke, K. A., Ces, O., Conn, C. E., Seddon, J. M. & Templer, R. H. (2007). Langmuir, 23, 7276–7285. Web of Science CrossRef PubMed CAS Google Scholar
This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.