SPIRE, Surface Projection Image Recognition Environment for bicontinuous phases: application for plastid cubic membranes

Bicontinuous membranes in cell organelles epitomise nature’s ability to create complex functional nanostructures. Like their synthetic counterparts, these membranes are characterised by continuous membrane sheets draped onto topologically complex saddle-shaped surfaces with a periodic network-like structure. In cell organelles, their structure sizes around 50–500 nm and fluid nature make Transmission Electron Microscopy (TEM) the analysis method of choice to decipher nanostructural features. Here we present a tool to identify bicontinuous structures from TEM sections by comparison to mathematical “nodal surface” models, including the hexagonal lonsdaleite geometry. Our approach, following pioneering work by Deng and Mieczkowski (1998), creates synthetic TEM images of known bicontinuous geometries for interactive structure identification. We apply the method to the inner membrane network in plant cell chloroplast precursors and achieve a robust identification of the bicontinuous diamond surface as the dominant geometry in several plant species. This represents an important step in understanding their as yet elusive structure-function relationship.


32
Biological membranes, dynamic yet stable unique assemblies of lipids and proteins, are selective 33 barriers and enzymatically active regions playing a crucial role in orchestrated cells' functioning. 34 From a structural point of view, they are described mainly as flat sheets or small folded isolated  plane. Note that the space between voxels is just for visualisation and is not existent in the simulation. (2) All voxels are evaluated using the mathematical model of the surface structure: the slice is superimposed in the oriented membrane structure; a voxel is marked if located within a membrane (2') The projection is computed by adding the values of voxels with identical and coordinates, i.e. all voxels congruent in viewing direction. Each marked voxel is valued as "1" , whereas unmarked voxels do not contribute (3) The simulated planar projection, a pixel image where each pixel brightness holds the number of marked, thus membrane, voxels. Its resolution is determined by the number of voxels ( , ) in the initial slice our tool for a large community of biologists. The intuitive SPIRE graphical user interface (GUI), see  Although we developed SPIRE to investigate cellular cubic membranes visualized in TEM images, 164 it is also a suitable tool for analyzing microscopy images of any highly symmetric arrangements 165 such as cuboids or polymer assemblies. In the latter, we see similar geometries to those found in  176 This article first describes the main concepts and algorithms of the software, followed by a 177 detailed walkthrough of a structure identification process using the PLB arrangements in etiolated 178 seedlings as an example. The Appendix contains a video tutorial presenting the tool and workflow, 179 as well as information about the built-in structures and parameters used throughout this article.

Numerical procedures
182 Figure 2 illustrates the major steps in the simulation process used in SPIRE. This software simulates 183 TEM images of ultrathin sections of biological samples, which are essentially planar projections of 184 the structure inside the sample. In general, a planar projection is a flat, thus 2D image of a 3D 185 structure. Though this transformation to a lower dimensional space causes a loss of information, 186 in many cases it yields an image providing a better overview than the 3D structure.

187
Here, we aim for a projection emulating a TEM image of a 3D structure. The latter encodes 188 the amount of beam attenuating material along its path that generates the microscope image: 189 whereas dark regions in the image represent areas of a large amount of material along the beam 190 path, brighter regions present less material.

191
To simulate those projections, a model of the sample is created by using one or multiple mem-192 branes, modeled as minimal or negatively curved surfaces. A discrete grid of points, where each 193 point is either marked as "attenuating", or not marked, that is "translucent", is then used to com-194 pute the planar projection. In this section details on the underlying processes and methods will (see Crystallographic nature of highly symmetric membranes for more details). Note, that due 216 to their repeating nature, all representations of triply-periodic surfaces can be expressed solely 217 using periodic trigonometric functions. Although the presented formula yields a surface with a 218 topology and geometry equivalent to the lonsdaleite structure, it is not a minimal surface. For the 219 implementation in the software we therefore used a numerically optimised version 1 . 220 1 A triangulated, minimal lonsdaleite surface was created using an input file by Ken Brakke for the surface evolver (Brakke, 1992(Brakke, , 2008. A voxelized version of the surface was then created by marking points on a rectangular grid as "1" on one side of the surface and "-1" on the other side. When interpolating, this grid has then the value "0" everywhere on the surface. This discrete, 3D test function was then approximated by a Fourier series ( , , ), using a Fast-Fourier-Transform (FFT). This 6 of 36 All TPMSs divide space into two intertwining channels. For this software, the nodal representa-   So far, only systems related to the double "gyroid" (or other "double phases") where a mem-

263
Fourier series then has a root everywhere on the surface and for ( , , ) is thus a nodal representation of the surface. The same numerical protocol was used to compute nodal representations of two cubic rod packings, namely the -Mn and Σ + rod packings. The original input was generated by placing cylinders of a given radius along the invariant axes of the rods, as described in O'Keeffe et al. (2001). Note that for convenience, in this implementation, instead of using the canonical choice of lattice vectors for a hexagonal symmetry, we use orthogonal lattice vectors with a rhomboidal symmetry. Please refer to the appendix Fundamental unit cells of built-in surfaces for the detailled information. 2 Whereas the surface ( , , ) = 0 will converge towards the true minimal surface by adding more terms to the series expansion, this is not the case for surfaces with ( , , ) = where ≠ 0. The latter, however are topologically equivalent (within a symmetric interval of values) non-minimal, triply periodic surfaces.

Figure 3. Multi-layer membrane structures and channel enumeration
A primitive surface multi-layer membrane system and its 2D projection in (100) orientation with two membranes of width at a distance of computed as parallel surfaces from the level-set membrane, the minimal surface at ( , , ) = 0, shown in red. The latter is only computed internally and does not show in the projection. The inner membrane is inside of the level-set membrane, therefore has a negative distance. The numbers denote the channel numbers of a total of 5 channels, of which 3 are "true" channels and 2 are membranes, also considered as channels internally.
brane (modelled as absorbing) separates two aqueous channels (modelled as translucent) were  Analogously, the software allows for membranes to be marked translucent to allow for more com- and will be filled with a discretized version of the membrane structure. To store the latter, the sim-281 ulation box is fitted with a regular, rectangular grid of a total of ⋅ ⋅ sites, see part 1 in Figure 2. 282 On each site a small cuboid with dimensions ( , , ), a so-called voxel, is placed, such as there 283 is no overlap or space between two adjacent voxels. Each of these voxels can be "unmarked", i.e. 284 translucent, or "marked", i.e. opaque. In Figure 2 this property is represented by the color "white" 285 and "red" and will be assigned in a later step. the resolution of the TEM image. 294 The current version of the software is designed to model to periodic membrane structures. The  The entire information of a periodic structure is thus stored in only a single translational UC 303 and the three replication directions, called lattice vectors , and . We conveniently choose these 304 vectors as the bounding edges of the UC (see Figure 4), although other choices exist and may be 305 more fitting for different purposes. As a consequence the size and shape of the latter is solely 306 determined by a choice of the three lattice vectors , and , i.e. the replication directions. defining the viewing angle on the structure as vector ℎ , as discussed above. However, since the 352 viewing angle in the tool is fixed to the axis, the vector n ℎ needs to be aligned with the z axis, 353 applying the same rotations on the membrane structure, as visualized in Figure 6. 354 The alignment is performed in three separate rotations and a final translation, as shown in where is the smallest integer, for which (ℎ ) is a triplet of integers without a common divisor. The polar angle Φ is the angle between the axis and the projection of ℎ onto the -plane, the azimuthal angle Θ is the angle between the axis and ℎ . Figure 6. Orienting the membrane structure in the virtual sample To simulate a projection with the desired viewing angle, the structure needs to be inside the slice. Starting at the orientation of the fundamental unit cell (UC), the structure is first rotated around the axis by the polar angle Φ, rotating the vector into the -plane, and subsequently around the axis by the azimuthal angle Θ, aligning with the axis and thus the viewing direction. The last rotation around the axis aligns the in-plane vector v with the x axis, such that the inclination UC has minimal volume given the normal vector n ℎ . In a last step the structure can be translated along the axis to choose the termination. ure 4 shows an example vector , which is long compared to the size of the fundamental UC. This 370 behavior directly translates to three dimensions.

371
Having the orientation fixed, the last step is to fix the translational degrees of freedom. This 372 can be imagined as the oriented slice being a stencil, cutting a rectangular piece out of the infinite, 373 periodic membrane structure at different locations. Moving a stencil along a periodic structure 374 does not significantly change its content, given the stencil is at least the size of a inclination UC: 375 the slice will in most cases contain one or multiple copies of the UC. As most biological samples 376 will be ultra-thin slices, i.e. have a much larger base than thickness, this caveat is met in most 377 cases in the and direction of the sample, i.e. its width and height, but not the direction, the 378 viewing direction. Thus, whereas translating the slice through the membrane structure in and 379 direction does barely affect the projection, moving it along the direction will impact the projection  The planar projection of a membrane structure is a pixel image, thus an array of × pixels, 388 each with a value indicating the brightness of the pixel. In the biological sample, a dark pixel in a 389 TEM image indicates that much of the brightness of the incident beam has been lost due to a large 390 amount of attenuation with matter, i.e. much electron-dense material is located along the path of 391 the beam corresponding to that particular pixel.

392
This behavior is imitated in SPIRE: each voxel marked as being inside a membrane is assigned 393 a numerical value of "1", whereas all unmarked pixels, thus representing aqueous phases, are 394 assigned a value of "0". A value of "1" means that a membrane, therefore attenuating material, the connecting path to be cut and thus renders the network non-percolating.

427
This software uses percolation analysis to compute the percolation threshold of a channel, that 428 is the maximum diameter of a body (e.g., a sphere or molecule) which can move freely through a 429 channel of the entire structure without getting stuck at narrow passages (Mickel et al., 2008). This iii) the orientation of a structure in the TEM sample. We discuss how to use these parameters for 446 proper and efficient matching.

447
Since the size of the TEM sample is known, the first parameters which can easily be fixed to start analysis tools (e.g., ImageJ). 458 The next step is to make an initial (educated) guess for the structure type and orientation. A   Table 1). In specific cases, the PLB can adopt  (Figure 12). Note that a significant increase in volume proportion alone, with-504 out changes in UC size, causes a rise in the membrane area packed in a given volume (Figure 12). 505 Therefore, when the PLB size is maintained between different genotypes, such a balanced PLB 506 network might store significantly larger amounts of membrane components, including enzymatic 507 proteins and galactolipids crucial for efficient etioplast-chloroplast transition.

508
In Figure 13 we present the exemplary utilization of the percolation threshold function.The  Table 1 17 of 36   . This variability also influences the structure's membrane packaging potential; balanced PLBs can accumulate more membrane components in the given volume; all parameters used to generate projections are listed in Table 1 Figure 13. Penetrability of the three-dimensional prolamellar network -combining ultrastructural and molecular data. To estimate the maximum diameter of a sphere that can freely penetrate through a channel of the entire structure, the "percolation threshold" implemented in SPIRE can be computed (A). Here we used an example of plastid ribosome (B -PDB 5MMM, spinach Spinacia oleracea 70S chloroplast ribosome modeled using Chimera software) which size is comparable with diameters of the oat Avena sativa prolamellar body (PLB) channels. Using the percolation threshold function, it is possible to estimate whether a molecule whose spatial structure is already revealed could move freely through the channels of a particular cubic surface. The percolation threshold (31.57 nm) of a given network (unit cell size 80 nm, volume proportion 0.25) , identified by a structure's projection (D), can be calculated in SPIRE. The diamond network of oat PLB in (111) direction is presented with ribosomes in the same scale, region marked with white square presents superposition of computed projection and TEM image using multiply blend mode. Note that proper calculations are also provided for the surfaces recognized by projections in which the channel's maximal diameter (40.11 nm) is not visible. In this example, ribosome size is smaller than the oat PLB's percolation limit, which indicates that this molecule can move freely through the larger aqueous channel of the PLB network (stroma) of oat, hypothetically fulfilling its biological function directly inside the cubic structure; all parameters used to generate projection are listed in Table 1.

Discussion
In this work, we introduced a new software tool (SPIRE), which, based on "nodal surface" mod-   With SPIRE's source code being designed to be usable as a library, as well as being available 564 freely under an open source license, we also see the possibility to integrate the code into existing 565 tools and software. This way, especially once the identification process has been fully automated, 566 we hope to make the identification of structures from stereographic images available to a large 567 22 of 36 and broad audience across many fields.

568
For now, simultaneously with this work, we provide a commonly available video tutorial (http:// 569 chloroplast.pl/spire) that provides a comprehensive guide for the efficient semi-automatic matching 570 procedure and presents all basic functionalities of SPIRE, the measurement tab in particular.

571
An anticipated broad interest in periodic structure recognition in biological samples enabled 572 by SPIRE will lead to the proper identification of many naturally occurring bicontinuous structures,  The code is designed to allow an easy implementation of further surface types, given in the   (1 1 1) 5