Abstract
Complementing collections of 3D time-lapse image data with comprehensive manual annotations is an extremely laborious and often impracticable task, which hinders objective benchmarking of bioimage analysis workflows as well as training of widespread deep-learning-based approaches. In this paper, we present a novel simulation system capable of generating synthetic 3D time-lapse sequences of multiple mutually interacting cells with filopodial protrusions, accompanied by inherently generated reference annotations, in order to stimulate the development of fully 3D bioimage analysis workflows for filopodium segmentation and tracking in complex scenarios with multiple mutually interacting cells. The system integrates its predecessor, which was designed for single-cell, collision-unaware scenarios only, with proactive, mechanics-based handling of collisions between multiple filopodia, multiple cell bodies, or their combinations. We demonstrate its potential on two generated 3D time-lapse sequences of multiple lung cancer cells with curvilinear filopodia, which visually resemble confocal fluorescence microscopy image data.
You have full access to this open access chapter, Download conference paper PDF
Similar content being viewed by others
Keywords
1 Introduction
The imaging of cell migration and inter-cellular interactions is essential to understand various physiological and pathological processes, such as tissue morphogenesis, wound healing, embryonic development, and cancer extravasation [2, 8, 9]. A remarkable role in these processes involving cell interactions is attributed to filopodia. These are cell membrane protrusions driven by dynamic bundles of aligned actin filaments, with not fully revealed detailed biology [10, 25].
The modern optical microscopy and fluorescent reporters facilitate observations of filopodial formation and dynamics in diverse three-dimensional microenvironments at unprecedented spatio-temporal resolutions [12]. However, the lack of robust and fully automatic bioimage analysis workflows for the quantification of filopodium-mediated processes in 3D+t image data impels biologists to perform such analyses solely on 2D+t image data in daily practice [1, 10, 22, 25].
A fundamental prerequisite to bridge this gap is the availability of diverse 3D+t image datasets, accompanied by reference annotations. Such datasets allow the bioimage analysis community to objectively benchmark developed bioimage analysis workflows. Nonetheless, they are difficult to obtain due to known subjectivity, high error-proneness, and extreme labor needed for manually annotating them [6]. Therefore, the focus has recently been put on developing simulation systems that routinely generate synthetic image data, accompanied by inherently generated reference annotations, that qualitatively and quantitatively resembles real counterparts [24]. Indeed, the synthetic image data was extensively used for benchmarking of diverse time-lapse bioimage analysis algorithms [4, 23] and for training of deep-learning approaches for various bioimage analysis tasks [3, 26].
In this paper, we build on a recent simulation system capable of generating 3D time-lapse sequences of single motile cells with filopodial protrusions of user-controlled structural and temporal attributes [19]. Namely, we extend this single-cell, collision-unaware system to be able to generate 3D time-lapse sequences of multiple mutually interacting motile cells. By introducing proactive, mechanics-based handling of collisions, the new system can tackle various types of cell-to-cell interactions, such as the contacts between multiple filopodia, multiple cell bodies, and their mutual combinations. Its capabilities are demonstrated by generating two synthetic 3D time-lapse sequences of two and three mutually interacting lung cancer cells with single-branch filopodia of fixed geometries, visually resembling real lung cancer cells acquired using a confocal fluorescence microscope.
2 Methods of the Simulation System
In this section, we describe our simulation system, including its overview and the methodology used for modeling deformations of multiple cells with filopodia, handling their collisions, and generating realistic, time-coherent cell textures.
2.1 System Overview
We consider a simulation of a set of cells where each cell is composed of an elastic cell body and a set of deformable filopodia of fixed geometries attached to the cell body surface. The initial geometries of the cell bodies and filopodia are generated using the geometry module of FiloGen [19]. As the cells move, collisions between them can naturally occur. The proposed simulation system is capable of handling collisions involving any pair of the cell components, such as two filopodia, two cell bodies, or a cell body and a filopodium, including auto-collisions which can potentially occur between two components of the same cell. Finally, the deformed cell geometries are exported in each step of the simulation and provided to the module generating time-coherent cell textures, as described in Sect. 2.4.
From the mechanical point of view, the simulation is based on constrained dynamics of elastic objects numerically implemented using a finite element method, with two types of constraints being modeled: equality constraints keep filopodia attached to the respective cell body surfaces, and inequality constraints guarantee non-penetration of cell components during collisions.
Each time step of the mechanical simulation involves three sequential phases: unconstrained dynamics, collision detection, and collision response. Initially, the motion of each cell component is modeled independently according to the Newton’s law of physics, which can result in violations of both equality and inequality constraints. Therefore, the vectors of constraint violations are computed by a proximity-based collision detection. Finally, the mechanics of each cell component is used to compute correction forces from the vector of violations and these forces are applied to the cell components to reimpose the constraints.
The entire simulated system can be regarded as an efficient domain decomposition where mechanical components of a large system are coupled using Lagrange multipliers computed with the Schur complement [5, 14].
2.2 Finite Element Models of Cell Components
The dynamics of each cell component c of the simulated system is given as
where \(\mathbf {q}_c, \dot{\mathbf {q}}_c, \ddot{\mathbf {q}}_c\) are respectively position, velocity, and acceleration of the mesh nodes representing the deformable cell component, \(\mathbf {M}_c\) is an inertia matrix, \(\mathbf {f}_c\) is a vector of external forces, and \(\mathbf {g}_c\) is a vector of internal elastic forces. Finally, the term \(\mathcal {H}_c^\top \varvec{\lambda }\) represents the mechanical constraints, with \(\varvec{\lambda }\) being correction forces derived from the collision response. The vector \(\varvec{\lambda }\) is defined in a constraint space shared by all colliding and coupled components. Therefore, it must be mapped to the component space through a generally nonlinear mapping \(\mathcal {H}_c\), represented by a rectangular matrix \(\mathbf {H}_c\) for simplicity. The internal forces are represented as a generally nonlinear function defined as
Its derivatives, \(\mathbf {K}_c=\frac{\partial \mathbf {g}_c}{\partial \mathbf {q}}\) and \(\mathbf {B}_c=\frac{\partial \mathbf {g}_c}{\partial \dot{\mathbf {q}}}\), are known as the stiffness and damping matrices, respectively.
The system described by (1) is integrated using the Euler implicit scheme: denoting t the previous time step, the velocity update in the actual time step \(t+h\) is computed as
being discretized using finite element models that reflect the cell component type. Linear co-rotational beam elements based on the Tymoshenko formulation [17] are considered in the case of a filopodium \(\varphi \). Organized in piecewise linear segments, its nodes \(\mathbf {q}_\varphi \) have six degrees of freedom (DoFs). The stiffness matrix \(\mathbf {K}_\varphi \), consisted of 12\(\times \)12 blocks as each beam relates two 6-DoF nodes, is parametrized with Young’s modulus \(E_\varphi \), Poisson ratio \(\nu _\varphi \), and the beam radius \(r_\varphi \). As for a cell body \(\beta \), linear tetrahedral P1 elements are considered and co-rotational formulation of linear elasticity is used [11]. The nodes \(\mathbf {q}_\beta \) have three DoFs and the stiffness matrix \(\mathbf {K}_\beta \), consisted of 12\(\times \)12 blocks relating the contributions of four 3-DoF nodes connected to a P1 element, is parametrized with Young’s modulus \(E_\beta \) and Poisson ratio \(\nu _\beta \). Finally, for both component types, the Rayleigh approximation of the damping matrix is employed: \(\mathbf {B}_c = r_K\mathbf {K}_c + r_M\mathbf {M}_c\), with \(r_K\) and \(r_M\) being the Rayleigh stiffness and the Rayleigh mass, respectively.
2.3 Dynamic Simulation with Collisions
Initially at \(t=0\), all simulated cells are in a collision-free configuration. Each filopodium \(\varphi \) is associated to a root \(\mathbf {r}_{\varphi \beta }\) which is located on the surface of the cell body \(\beta \) and initially, the distance \(|\mathbf {r}_{\varphi \beta }, \mathbf {q}_{\varphi 0}| = 0\) where \(\mathbf {q}_{\varphi 0}\) is the first node of the filopodium \(\varphi \). For given cell body \(\beta \), all the filopodium roots \(\mathbf {r}_{\varphi \beta }\) are geometrically mapped to the mesh of parent body \(\beta \): this mapping is implemented via barycentric coordinates of the root with respect to its embedding tetrahedron. The coordinates are computed at the initialization and remain constant during the simulation as described in [15]. Therefore, as the cell body moves, the filopodium roots strictly follow this motion.
Two sources of motions are applied to each cell: the center of each cell body is moved along a predefined trajectory, and the forces of predefined directions and amplitudes are applied to the filopodial tips. Then, given the actual configuration obtained in time t, the three phases of step \(t+h\) consist in the following:
Unconstrained Dynamics. For each component c, (3) is assembled using the last known configuration given by \(\mathbf {q}^t_c\) and \(\dot{\mathbf {q}}_c^t\), and the corresponding system is solved with the constraint forces \(\varvec{\lambda }= 0\). The resulting velocity corrections \(d\dot{\mathbf {q}}_c^\text {uncon}\) are integrated over the time step to get the unconstrained configuration given by the vector of positions \(\mathbf {q}_c^\text {uncon}\) of each component of the simulated system. The assembled systems are solved directly, using a highly optimized solver based on the Thomas algorithm for block tridiagonal matrices [21] in the case of filopodia, and using the Pardiso solver [16] in the case of cell bodies. These solvers also construct the decomposition of the system matrix \(\mathbf {A}_c\) needed to compute the Schur complement \(\mathbf {H}_c \mathbf {A}_c^{-1} \mathbf {H}_c^\top \) later.
Collision Detection. During the unconstrained motion, individual cell bodies and their filopodia follow the Newton’s law of physics as independent entities, which can possibly lead to the violations of equality and inequality constraints.
In the case of equality constraints, for each filopodium \(\varphi \) and its parent body \(\beta \), the first node of the filopodium \(\mathbf {q}_{\varphi 0}^\text {uncon}\) does not have to necessarily coincide with the current position of the corresponding filopodium root \(\mathbf {r}_{\varphi \beta }^\text {uncon}\). Therefore, the violation is defined as the difference \(\delta _E^\text {uncon} = |\mathbf {r}_{\varphi \beta }^\text {uncon} - \mathbf {q}_{\varphi 0}^\text {uncon}|\), with E running over all pairs of \(\varphi \) and \(\beta \) present in the simulation.
The violations of inequality constraints are computed by a proximity-based collision detection directly available in the SOFA framework [7]. It produces a vector of violations \(\delta _I^\text {uncon}\), with I running over all inter-penetrations between collision primitives associated to each filopodium and each cell body. The vector values correspond to the negative depths of individual inter-penetrations along the axis and surface normals for filopodia and cell bodies, respectively.
Besides the violations of both types of constraints, the corresponding rectangular matrices \(\mathbf {H}_c\) are constructed for each component c, following the protocol introduced in [5]. In the case of equality constraints, the matrix contains the barycentric coordinates of the filopodium roots, and it maps the exact contact positions as they have occurred between particular collision primitives in the case of inequality constraints.
Collision Response. The vector of constraint forces \(\varvec{\lambda }\) (one force for each equality and inequality constraint) is to be found by solving (4). In the case of equality constraints, the constraint is satisfied if \(\delta _E = 0\). In the case of inequality constraints representing frictionless contacts, the constraint is satisfied when it respects Signorini’s law given as \(\delta _I \le 0 \perp \varvec{\lambda }_I \ge 0\), representing the complementarity of contact forces \(\varvec{\lambda }_I\) and positive gap \(\delta _I\) between two objects. In other words, the inequality constraint is satisfied when either the negative violation becomes a positive gap (i.e., no contact occurs, so \(\varvec{\lambda }_I = 0\)) or non-zero contact forces are applied to remove the inter-penetration and then \(\delta _I = 0\). The relation between the violations and correction forces is based on the Schur complement:
where c and d denote two components of the simulated system which are coupled by given constraint, and \(*\) denotes either E or I subscripts. The c-d pair is always formed by one filopodium and one cell body for equality constraints, and any pair of colliding objects for inequality constraints.
Putting together the equations from all equality and inequality constraints, we arrive at a linear complementarity problem. We solve for \(\varvec{\lambda }\) using the projected Gauss-Seidel method [5]. Finally, the found vector of correction forces is used to perform the final set of corrections of each component, \(\mathbf {q}_c^{t+h} = \mathbf {q}_c^\text {uncon} + \mathbf {A}_c^{-1} \mathbf {H}_c \varvec{\lambda }\).
2.4 Generation of Cell Texture
Unlike the simulations limited to 2D image data, where the developers of simulation platforms have to cope with possibly overlapping cells, we model the cells fully in 3D space. These cells can touch one another, but they cannot mutually penetrate. Therefore, the texture of each cell is processed separately, following the protocol introduced in [18, 19]. The texture of each simulated cell is represented as a cloud of many simulated fluorescent particles. They form a phantom image in which all cells are presented and in which the intensity of every voxel is proportional to the number of fluorescent particles that occupy the voxel volume. The initial distribution of these particles is taken from an image with an isotropic mosaic of procedural texture [13], being restricted to the cell volume. When rendering successive frames of the time-lapse sequence, the particle positions follow the cell geometry changes by tracking displacements of the specific mesh vertices. To obtain the final image that looks as if acquired using a fluorescence microscope, the phantom image is submitted to a virtual microscope [20]. In particular, the phantom image is blurred by an experimental point spread function of the real optical system, unevenly illuminated, resampled to the resolution of the simulated microscope setup, and finally degraded by photon shot noise, dark current, and read-out noise.
3 Results
To demonstrate the capability of our simulation system, we imitated the interactions of CRMP-2-overexpressing lung cancer cells, being transiently transfected with the green fluorescent protein conjugated to actin and exhibiting mesenchymal migration typical for cancer cell extravasation [9]. As a reference, we considered real 3D+t image data with the voxel size of \(0.126\times 0.126\times 1.0\,{\upmu }\hbox {m}\), being acquired every two minutes using a spinning disk confocal microscope equipped with a water Plan-Apo \(63\times /1.20\) objective lens. The cells had variable cell body textures and protruded very short, single-branch filopodia of low motility.
We performed two simulations consisting of two and three mutually interacting motile cells, respectively, applying in each the time step of \(h = 0.01\,\hbox {s}\), Young’s moduli of \(E_\beta = 10\,\hbox {kPa}\) and \(E_\varphi = 500\,\hbox {kPa}\), Poisson ratios of \(\nu _\beta = 0.45\) and \(\nu _\varphi = 0.4\), the Rayleigh mass and stiffness of 0.1 for both the cell bodies as well as filopodia, and the tolerance of \(10^{-6}\) for the constraint solver used in (4). The shift of each cell body centroid over one time step was within \(0.03\,{\upmu }\hbox {m}\), and the forces applied to the filopodial tips ranged from 100 to \(200\,{\upmu }\hbox {N}\). The resulting 3D time-lapse sequences are shown in Fig. 1, demonstrating time-coherent, mechanics-based mutual interactions between the cell bodies, filopodia, and also between the cell bodies and the filopodia. In order to qualitatively compare the synthetic cell textures with its real counterparts, an example of four successive XY and XZ slices of both is depicted in Fig. 2. Their quantitative comparison was already given in [19], together with a more detailed description of the protocol that we followed here when generating the textures of multiple cells. The two complete synthetic sequences, accompanied by the inherently generated ground truth, and the source codes, developed to generate these two sequences, are made publicly available, free of charge for any research and noncommercial purposes, at http://cbia.fi.muni.cz/research/simulations/multiple-cells-filopodia.html.
An example of four successive XY and XZ slices of the real (top) and synthetic (bottom) image data of two colliding cells. The voxel size is \(0.126\times 0.126\times 1.0\,{\upmu }\hbox {m}\). The white ticks show the positions of the other two orthogonal cross-sections of image data. The scale bars correspond to \(5\,{\upmu }\hbox {m}\).
In order to provide the reader with an estimate of how long the generation of the synthetic image data took, we measured the execution time of both simulations on a common workstation (Intel Xeon QuadCore 2.83 GHz, 32 GB RAM, Ubuntu 16.04 LTS). On average, the generation of one initial cell geometry took less than a minute. The modeling of cell motion and deformation, including handling of cell interactions, took roughly one and three minutes for the scenarios involving two and three cells, respectively. The image synthesis of 71 frames long time-lapse sequences of two and three cells, with the frame size of \(336\times 275\times 18\) and \(284\times 375\times 18\) voxels, respectively, took 56 and 84 min, respectively.
4 Conclusion
We have presented a simulation system that can generate synthetic 3D time-lapse sequences of multiple mutually interacting motile cells with filopodial protrusions of fixed geometries, accompanied by reference annotations. The system extends its predecessor [18, 19], designed for single-cell, collision-unaware scenarios solely and proven to generate realistic cell textures, by adding proactive, mechanics-based handling of collisions to naturally deal with the interactions between multiple cell bodies and filopodia. Along with the cell kinetics and material parameters provided, the system allowed us to generate two plausible-looking synthetic 3D time-lapse sequences of multiple lung cancer cells with single-branch filopodia. Showing time-coherent, mechanics-constrained development, and being accompanied by inherently generated ground truth, such synthetic 3D+t image data can be exploited for objective benchmarking of bioimage analysis workflows as well as for training of deep-learning approaches that target fully 3D filopodium segmentation and tracking for multiple mutually interacting cells. In future work, we intend to extend our system with bleb-like protrusions and filopodia dynamically evolving their geometries over time, along with investigating an appropriate dynamic model for entire cell kinetics and its subsequent quantitative analysis.
References
Barry, D.J., Durkin, C.H., Abella, J.V., Way, M.: Open source software for quantification of cell migration, protrusions, and fluorescence intensities. J. Cell Biol. 209(1), 163–180 (2015)
Biswas, K.H., Zaidel-Bar, R.: Early events in the assembly of E-cadherin adhesions. Exp. Cell Res. 1(358), 14–19 (2017)
Castilla, C., Maška, M., Sorokin, D.V., Meijering, E., Ortiz-de-Solorzano, C.: Segmentation of actin-stained 3D fluorescent cells with filopodial protrusions using convolutional neural networks. In: IEEE International Symposium on Biomedical Imaging, pp. 413–417 (2018)
Chenouard, N., et al.: Objective comparison of particle tracking methods. Nat. Methods 11(3), 281–289 (2014)
Courtecuisse, H., Allard, J., Kerfriden, P., Bordas, S.P.A., Cotin, S., Duriez, C.: Real-time simulation of contact and cutting of heterogeneous soft-tissues. Med. Image Anal. 18(2), 394–410 (2014)
Coutu, D.L., Schroeder, T.: Probing cellular processes by long-term live imaging-historic problems and current solutions. J. Cell Sci. 126(17), 3805–3815 (2013)
Faure, F., et al.: SOFA: a multi-model framework for interactive physical simulation. In: Payan, Y. (ed.) Soft Tissue Biomechanical Modeling for Computer Assisted Surgery. SMTEB, vol. 11, pp. 283–321. Springer, Berlin (2012). https://doi.org/10.1007/8415_2012_125
Haeger, A., Wolf, K., Zegers, M.M., Friedl, P.: Collective cell migration: guidance principles and hierarchies. Trends Cell Biol. 25(9), 556–566 (2015)
Jacquemet, G., Hamidi, H., Ivaska, J.: Filopodia in cell adhesion, 3D migration and cancer cell invasion. Curr. Opin. Cell Biol. 31(10), 23–31 (2015)
Jacquemet, G., et al.: FiloQuant reveals increased filopodia density during breast cancer progression. J. Cell Biol. 216(10), 3387–3403 (2017)
Nesme, M., Payan, Y., Faure, F.: Efficient, physically plausible finite elements. In: Eurographics, pp. 77–80 (2005)
Ortiz-de-Solórzano, C., Muñoz-Barrutia, A., Meijering, E., Kozubek, M.: Toward a morphodynamic model of the cell. Signal Process. Mag. 32(1), 20–29 (2015)
Perlin, K.: An image synthesizer. In: SIGGRAPH, pp. 287–296 (1985)
Peterlík, I., et al.: Fast elastic registration of soft tissues under large deformations. Med. Image Anal. 45(4), 24–40 (2018)
Peterlík, I., Duriez, C., Cotin, S.: Modeling and real-time simulation of a vascularized liver tissue. In: Ayache, N., Delingette, H., Golland, P., Mori, K. (eds.) MICCAI 2012. LNCS, vol. 7510, pp. 50–57. Springer, Heidelberg (2012). https://doi.org/10.1007/978-3-642-33415-3_7
Petra, C.G., Schenk, O., Lubin, M., Gärtner, K.: An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization. SIAM J. Sci. Comput. 36(2), C139–C162 (2014)
Przemieniecki, J.S.: Matrix structural analysis of substructures. Am. Inst. Aeronaut. Astronaut. J. 1(1), 138–147 (1963)
Sorokin, D.V., Peterlík, I., Ulman, V., Svoboda, D., Maška, M.: Model-based generation of synthetic 3D time-lapse sequences of motile cells with growing filopodia. In: IEEE International Symposium on Biomedical Imaging, pp. 822–826 (2017)
Sorokin, D.V., et al.: FiloGen: a model-based generator of synthetic 3D time-lapse sequences of single motile cells with growing and branching filopodia. IEEE Trans. Med. Imaging (2018, in press). https://doi.org/10.1109/TMI.2018.2845884
Svoboda, D., Kozubek, M., Stejskal, S.: Generation of digital phantoms of cell nuclei and simulation of image formation in 3D image cytometry. Cytom. Part A 75A(6), 494–509 (2009)
Thomas, L.H.: Elliptic problems in linear difference equations over a network. Watson Science Computer Lab Report, Columbia University, New York (1949)
Tsygankov, D., Bilancia, C.G., Vitriol, E.A., Hahn, K.M., Peifer, M., Elston, T.C.: CellGeo: a computational platform for the analysis of shape changes in cells with complex geometries. J. Cell Biol. 204(3), 443–460 (2014)
Ulman, V., et al.: An objective comparison of cell-tracking algorithms. Nat. Methods 14(12), 1141–1152 (2017)
Ulman, V., Svoboda, D., Nykter, M., Kozubek, M., Ruusuvuori, P.: Virtual cell imaging: a review on simulation methods employed in image cytometry. Cytom. Part A 89(12), 1057–1072 (2016)
Urbančič, V., et al.: Filopodyan: an open-source pipeline for the analysis of filopodia. J. Cell Biol. 216(10), 3405–3422 (2017)
Yao, Y., Smal, I., Meijering, E.: Deep neural networks for data association in particle tracking. In: IEEE International Symposium on Biomedical Imaging, pp. 458–461 (2018)
Acknowledgements
This work was financially supported by the Czech Science Foundation [GJ16-03909Y to I.P., D.V.S., and M.M.], the Russian Science Foundation [17-11-01279 to D.V.S.], and the German Federal Ministry of Research and Education (BMBF) [031L0102 (de.NBI) to V.U.].
Author information
Authors and Affiliations
Corresponding author
Editor information
Editors and Affiliations
Rights and permissions
Copyright information
© 2018 Springer Nature Switzerland AG
About this paper
Cite this paper
Peterlík, I., Svoboda, D., Ulman, V., Sorokin, D.V., Maška, M. (2018). Model-Based Generation of Synthetic 3D Time-Lapse Sequences of Multiple Mutually Interacting Motile Cells with Filopodia. In: Gooya, A., Goksel, O., Oguz, I., Burgos, N. (eds) Simulation and Synthesis in Medical Imaging. SASHIMI 2018. Lecture Notes in Computer Science(), vol 11037. Springer, Cham. https://doi.org/10.1007/978-3-030-00536-8_8
Download citation
DOI: https://doi.org/10.1007/978-3-030-00536-8_8
Published:
Publisher Name: Springer, Cham
Print ISBN: 978-3-030-00535-1
Online ISBN: 978-3-030-00536-8
eBook Packages: Computer ScienceComputer Science (R0)