3D immuno-confocal image reconstruction of fibroblast cytoskeleton and nucleus architecture

Computational models of cellular structures generally rely on simplifying approximations and assumptions that limit biological accuracy. This study presents a comprehensive image processing pipeline for creating unified three-dimensional (3D) reconstructions of the cell cytoskeletal networks and nuclei. Confocal image stacks of these cellular structures were reconstructed to 3D isosurfaces (Imaris), then tessellations were simplified to reduce the number of elements in initial meshes by applying quadric edge collapse decimation with preserved topology boundaries (MeshLab). Geometries were remeshed to ensure uniformity (Instant Meshes) and the resulting 3D meshes exported (ABAQUS) for downstream application. The protocol has been applied successfully to fibroblast cytoskeletal reorganisation in the scleral connective tissue of the eye, under mechanical load that mimics internal eye pressure. While the method herein is specifically employed to reconstruct immunofluorescent confocal imaging data, it is also more widely applicable to other biological imaging modalities where accurate 3D cell structures are required.


| INTRODUCTION
In vivo, cells transmit mechanical signals into intracellular responses via rearrangement of the cytoskeleton. Understanding the precise contribution of the cytoskeleton to cell biomechanics is vital to model cell and tissue-level strains that control biological homeostasis and its perturbation in disease. The complexity and dynamics of the cytoskeleton have proven to be the main hurdle in developing numerical models, with myriad oversimplifications and exclusion of key components being common. Current cell models make a variety of simplifications, from depicting the cell as an ellipsoid to containing only random fibres. [1][2][3] Whilst several earlier publications have utilised cell-specific meshes for computational analysis, they, however, did not accommodate the proper cytoskeletal orientation and commonly depict only a singular cytoskeletal network. [4][5][6][7] Modern-day microscopy systems allow production and recording of extensive volume image data sets though non-destructive optical sectioning. In combination with proprietary software packages, the captured information can be reconstructed in three-dimensional (3D) space and further exported and refined. Fibroblasts are the sole cellular contributors in regulating the mechanics of the peripapillary sclera in response to intraocular pressure (IOP), a process, which underpins the blinding disease glaucoma. 8 In this study, bovine scleral fibroblasts were subjected to an in vitro loading regime that mimicked physiological IOP strains, and cytoskeletal element organisation was visualised post-load using immunofluorescent labelling and laser scanning confocal microscopy. We detail an approach for generating accurate representations of the three cytoskeletal components (both fibre and filament networks) and nuclei from collected experimental image data sets which can be utilised in modelling the mechanical behaviour of cells. Furthermore, to the best of our knowledge this is the first work to create cell-specific digital reconstructions of the vimentin intermediate filaments and β-tubulin microtubular networks.

| MATERIALS AND METHODS
Reagents were obtained from Sigma-Aldrich (Poole, UK) unless specified otherwise.

| Fibroblast mechanical loading
To mimic physiological loading, the bovine eye intraocular pressure (IOP) was assigned 27 +/− 5.4 mmHg (3.6 +/− 0.72 kPa). 9 The Young-Laplace equation (σ = pr/2 t) was utilised to determine IOP-induced hoop stress (σ) after approximating eye globe shape as spherical with an average tissue thickness (t) of 1.6 mm and a radius (r) of 16.2 mm, and a mean Young's modulus (E) for bovine eyes of 9 MPa. 10 As strain (ε) is the ratio of stress divided by Young's modulus (ε = σ/E), the physiological tensile strain applied was between 0.26% and 1.8%.
Physiological cyclic tensile strain (CTS; 0.26%-1.8% strain) was globally applied to all cells at 37 C (FX3000 system, Flexcell International); unstrained fibroblasts served as controls. Cells were strained for 1 hour using a 1 Hz frequency, representative of the in vivo IOP pulse frequency. 11 Fibroblasts were fixed 1 hour after CTS ceased using 2% (w/ v) paraformaldehyde for 15 minutes. F-actin stress fibre anisotropy was quantified (mean ± standard deviation) using the ImageJ plugin FibrilTool, 12 with higher values (up to a maximum of 1) indicating a greater fibre alignment, whilst lower values (down to a minimum of 0) are closer to a completely random distribution.
Cytoskeletal structures were imaged using a Zeiss LSM880 Airyscan™ upright laser scanning confocal microscope operated through Zen Black software (Carl Zeiss Ltd, Germany). Image acquisition was performed using a Plan-Apochromat 63x/1.4 Oil DIC M27 oil immersion objective lens. Laser line switching was used for sequential scanning of DAPI (Ex max: 358 nm; Em max: 461 nm) and Alexa 488 (Ex max: 495 nm; Em max: 520 nm) in Airyscan Fast mode. The scan area was 134.95 μm by 134.95 μm, with a pixel size of 0.0706 × 0.0706 μm (1912 × 1912 pixels) and a pixel dwell time of 0.55 μs. Z-stacks of 8-bit optical sections (recorded in CZI Zeiss file format) were taken through the entire depth of individual fibroblasts using 2× Nyquist sampling resulting in a Z-step of 0.14 μm (average of 32 optical slices per Z-stack). 2D Z-stack montages were created in ImageJ/Fiji. 13 Representative cells were imaged across three different wells per plate (n = 3 wells) across three independent experiments (N = 3) accounting for a total of 129 to 135 cells per label per regime.

| 3D reconstruction of cytoskeletal components and nuclei
Image stacks of cytoskeletal networks and nuclei were reconstructed to 3D isosurfaces using Imaris 9.2 (Bitplane, Switzerland); background subtraction thresholding was enlisted, subtracting variable background intensity from each voxel, followed by automatic tracing and automatic surface rendering, based on an algorithm by Costes and Lockett. 14 Resultant 3D isosurfaces were exported as virtual reality modeling language file format (WRL) file format. Tessellation simplification, carried out in MeshLab (December 23, 2016 built), reduced the number of composing elements of the initial meshes. 15 Quadric edge collapse decimation (QECD) was applied to the entire geometries of imported WRL files, with selected parameters for preserving mesh topology boundaries, and automatic vertex removal and repositioning. To ensure preservation of the reconstructions the cytoskeletal network polygon reduction was set to 5% of the original number of elements, and 10% for the nuclei, on the basis of testing the best overall topology contour retention with lowest number of vertices.
Geometries were exported from MeshLab as OBJ files, remeshed and made uniform in Instant Meshes. 16 Element shape remained triangular to better suit the curved topology of cytoskeletal fibres and nuclei. Target vertex number was set as close as possible to the input geometry number of vertices. The OBJ files were re-loaded in MeshLab to check for topology inaccuracies (such as holes or semi-detached regions) and exported to STL (stereolithography) file format. Any mesh inconsistencies were resolved in Netfabb 7.4 (Netfabb GmbH, Germany) to ensure an intact mesh. Volume meshing was accomplished through either of two pathways: (I) STL files were imported to Gmsh 4.3.0 where volume was produced from continuum first-order C3D4 tetrahedron elements, composed of four vertices forming four triangles between them. 17 The 3D mesh was exported as INP, a model-associated format of the finite element solver ABAQUS 6.14/ complete abaqus environment (CAE) (Dassault Systèmes Simulia Corp). (II) STL files were directly imported to ABAQUS as an orphan mesh through the 3D Mesh to Geometry plugin, with subsequent volume rendering.

| RESULTS AND DISCUSSION
A pipeline to accurately generate cell-specific 3D modelling reconstructions of cytoskeletal components and nuclei is presented (Figure 1). First, the input confocal image Z-stack (Figure 2A) voxel position and intensity details are utilised in Imaris to reproduce an isosurface triangulation ( Figure 2B) on the basis of the Marching Cubes algorithm. [18][19][20] Initial geometrical data obtained from Imaris is oversampled, containing low-quality triangular elements with sharp angles. This is due to the software preserving the smooth traced topology of curved fibres without minimising the number of vertices produced, which can lead to processing issues later. Furthermore, the initial meshes are often uneven and if utilised directly for modelling may lead to unreliable results. 21 As the performance of modelling analyses, such as finite element (FE) approaches, is critically dependent on input mesh quality, the initial isosurfaces are unsuitable in this state and require optimisation procedures. Exporting Imaris files in MeshLab allows for merging of duplicate vertices or removing those with quality lower than an assigned threshold. 15 In instances of geometries composed of millions of elements, specific isosurface triangulation decimation permits reduced complexity; the output polygon mesh contains a smaller number of elements and vertices, whilst specific properties of the original mesh for example, geometric distance between vertices and visual appearance of the mesh, can be controlled by user-assigned criteria. QECD is such an algorithm, which also employs penalising parameters to control mesh topology retention ( Figure 2C). 22 Following decimation, sufficient numbers of elements need to be preserved as large reductions can impact on structural composition.
Poor element shape (high aspect ratio) can affect convergence of the mesh, causing significant errors in numerical simulations, 21 making it preferable for the triangulation to be composed of identical equilateral triangles. Therefore, to make the geometry uniform, a retopology, or remeshing, is required. Instant Meshes is used to produce an isotropic mesh of triangular (Supplementary Figure 1) or quadratic elements. 16 Remeshing also includes assignment of a target vertex count due to the uniformity of the refined mesh, allowing for only a certain number of final vertices. Triangular geometries are exported to Gmsh, from MeshLab, as STL files, and subsequently converted to INP files ( Figure 2D,E). 17 This is performed to add volume to the triangulations as they only possess surface information at that point. After mesh retopology and refinement, Gmsh can add an internal volume composed of tetrahedral elements. Tetrahedrons containing surface triangular elements will be affected by their heterogeneity, with large and uneven distributions causing many distorted tetrahedral elements, further necessitating the described simplification steps. An alternative approach for mesh volume generation is to directly import the refined STL file into ABAQUS, converting it to a geometry model using the 3D Mesh to Geometry ABAQUS plugin. The final workflow step involves loading the isotropic topology volume meshes into the FE solver ABAQUS.
To evaluate the pipeline, bovine scleral fibroblasts were selected as an in vitro model system to assess the effect of tensile strain on cytoskeletal architecture and nucleus morphology. Resultant confocal image data of scleral fibroblasts, subjected to physiological strain (0.26%-1.8%, based on normal IOP in the human eye) or unstrained, and immunolabelled for F-actin (microfilaments), vimentin (intermediate filaments) and β-tubulin (microtubules) were processed (Figure 1) with confocal image input and isotropic mesh output presented (Figure 3). Physiological strain resulted in both cell and nucleus elongation, concomitant with increased F-actin stress fibre alignment ( Figure 3B). Loaded fibroblasts exhibited a 225% increase in anisotropy (0.283 ± 0.076, n = 129 cells) compared to unloaded cells (0.126 ± 0.038, n = 135 cells; P < .001). These alterations, and the noticeable heterogeneous arrangement, were accurately reflected in the reconstructed meshes, demonstrating the reliability of the approach ( Figure 3D). The illustrated meshes better represent the cellular components, as they are constructed from experimental results rather than approximations. The triangulated geometries closely follow the contours and dimensions of the cytoskeletal network, with the selected rendering and meshing parameters and fibre density determining the overall reconstruction precision. This allows effective visualisation and rendering of individual fibres and fibre bundles, the latter not having previously been included in published models.
An obvious potential downstream application of the presented workflow is in FE modelling. Cytoskeletal mesh precision is key for cell-focussed FE approaches, as the cytoskeleton is integral for the cell's organisation, structure and biomechanical responses. 23 Although earlier publications utilised cell-specific meshes for FE analysis to represent the cellular shape and nucleus, they did F I G U R E 1 Workflow schematic for generation of triangular element geometries from image volume dataset. Software packages are highlighted by blue rectangles, whilst model steps are underlined. Green arrows indicate the main standard procedure pipeline, with the red arrow denoting an ABAQUS plugin dependent step not accommodate detailed cytoskeletal orientation. [4][5][6] The approach presented reduces the number of approximations by providing a reproducible workflow for 3D cellular reconstructions from imaging data. To the best of our knowledge this is the first study to create digital reconstructions of the other cytoskeletal networks that is, vimentin intermediate filaments and β-tubulin microtubule networks, in addition to previously reported F-actin F I G U R E 2 Representative images of a nucleus, F-Actin stress fibres, vimentin and β-tubulin cytoskeleton at different procedural steps. A, Maximum intensity projection of confocal image Z-stacks taken at 63× magnification and resolution of 1912 × 1,912 pixels; B, Surface triangulation in Imaris 9.2; C, Rendering in MeshLab; D, Initial mesh in ABAQUS 6.14/CAE; E, Simplified and uniform mesh in ABAQUS 6.14/CAE F I G U R E 3 Confocal microscopy image data of a nucleus, F-Actin stress fibres, vimentin and β-tubulin cytoskeleton (a, b) and corresponding reconstructed meshes (c, d) in ABAQUS 6.14/CAE. A, Confocal images with no applied mechanical load; B, Confocal images 1 hour after cessation of mechanical load (0.26%-1.8% strain, 1 Hz for 1 hour); C, Corresponding reconstructed meshes of Figure A; D, Corresponding reconstructed meshes of Figure B. Scale bar = 20 μm stress fibres. 7 Furthermore, for the first time, this study produces experiment-based geometries for all three cytoskeletal elements. As such, in representing an integrated model system of the cytoskeleton and the mechanically linked nuclei, it will be useful in further understanding the complex, coordinated dynamics among the individual cellular components. Moreover, the downstream applications of these reconstructions are not solely limited to numerical analyses such as FE approaches, but could also potentially be utilised in other avenues, such as 3D-bioprinting. 24,25 As the cytoskeleton is reconstructed as a single object, individual fibres are merged at their base making it difficult to assess whether the fibres were originally connected or positioned above each other, which may result from the voxel intensity rendering. Labelling for specific cytoskeletal accessory molecules for example, actin-binding proteins, would assist in determining their connectivity and positioning in denser fibre regions. Apart from vimentin intermediate filaments enveloping the nucleus (Supplementary Figure 2), the imaged fibres were mostly in the same plane as fibroblasts are relatively flat cells. A proper evaluation of other cytoskeletal parameters, such as fibre length, width and density would require higher resolution imaging methods beyond the scope of the current study.
The ability to study cellular and cytoskeletal deformations in response to applied loading is fundamental for understanding how deformations may damage cells and tissues and contribute to cell-driven tissue growth and remodelling. As the cytoskeleton is highly dynamic and in a constant state of reassembly, numerical simulations necessitate inclusion of actin stress fibre contractility and reshaping to produce realistic biomechanical responses. 1,26 A partial solution could be accomplished by utilising cytoskeletal meshes from different time points and mechanical loading conditions as references for dynamic cytoskeletal changes ( Figure 3). The Zeiss 880 confocal system utilised in this study comprises an Airyscan FAST module which is designed for fast super-resolved live cell imaging (resolution 1.7× that of conventional confocal optics and a 4× speed enhancement i.e., < 100fps). This system has previously been used to image dynamic changes in microtubular architecture in cardiac myocytes, 27 demonstrating its capability to capture and reconstruct sequences of confocal images to model cytoskeletal dynamics.
Future work envisions performing parallel immunoconfocal and second-harmonic generation imaging on loaded cells in a 3D culture environment to capture both the cellular and extracellular matrix (ECM) response, respectively, in order to up-scale the system to capture tissue-scale dynamics and model cell-ECM interactions in a closer to physiological situation. Such an approach has been successfully used to assess cellular mechanical response in other tissues, including cornea, 28 cartilage 29 and bone. 30

| CONCLUSIONS
The work here presents, to the best of our knowledge, the first experimentally determined cell-specific cytoskeletal models that include all major individual filament and fibre components. The refinement workflow presented can be applied to any imaging technique that produces volume data sets for example, confocal microscopy. The initial rendered geometries are controllably reduced, retaining topology whilst decreasing the number of elements, and then made isotropic, facilitating mesh quality for subsequent FE analysis. The results from two different mechanical loading conditions were acquired and reconstructed as volume meshes, which potentially could be implemented in inverse FE modelling to analyse cytoskeletal and nuclear mechanical properties. Such an approach can facilitate the study of diseases with biomechanical aspects, such as glaucoma, and has widespread application where modelling physiological and pathological conditions to identify disease progression is pivotal.