FiloGen: A Model-Based Generator of Synthetic 3D Time-Lapse Sequences of Single Motile Cells with Growing and Branching Filopodia Dmitry V. Sorokin, Igor Peterlík, Vladimír Ulman, David Svoboda, Tereza Nečasová, Katsiarina Morgaenko, Lívia Eiselleová, Lenka Tesařová, and Martin Maška Member, IEEE

Abstract-The existence of diverse image datasets accompanied by reference annotations is a crucial prerequisite for an objective benchmarking of bioimage analysis methods. Nevertheless, such a prerequisite is arduous to satisfy for time-lapse, multidimensional fluorescence microscopy image data, manual annotations of which are laborious and often impracticable. In this paper, we present a simulation system capable of generating 3D time-lapse sequences of single motile cells with filopodial protrusions of user-controlled structural and temporal attributes, such as the number, thickness, length, level of branching, and lifetime of filopodia, accompanied by inherently generated reference annotations. The proposed simulation system involves three globally synchronized modules, each being responsible for a separate task: the evolution of filopodia on a molecular level, linear elastic deformation of the entire cell with filopodia, and the synthesis of realistic, time-coherent cell texture. Its flexibility is demonstrated by generating multiple synthetic 3D time-lapse sequences of single lung cancer cells of two different phenotypes, qualitatively and quantitatively resembling their real counterparts acquired using a confocal fluorescence microscope.

I. INTRODUCTION

F ILOPODIA are plasma membrane protrusions governed by dynamic bundles of aligned actin filaments. Acting as guidance sensors and promoters of cell locomotion, they play crucial roles in many physiological and pathological processes, such as cell adhesion, chemotaxis, wound healing, embryonic development, or cancer extravasation [START_REF] Mattila | Filopodia: Molecular architecture and cellular functions[END_REF]- [START_REF] Jacquemet | Filopodia in cell adhesion, 3D migration and cancer cell invasion[END_REF].

Ever-growing technological advances in optical microscopy and fluorescent reporters have facilitated observations of these processes at unprecedented spatio-temporal resolutions [START_REF] Ortiz-De-Solórzano | Toward a morphodynamic model of the cell[END_REF], at the expense of increasing amounts of acquired image data to be processed and quantitatively evaluated, calling for robust and fully automatic bioimage analysis solutions [START_REF] Cardona | Current challenges in open-source bioimage informatics[END_REF], [START_REF] Carpenter | A call for bioimaging software usability[END_REF]. However, objective assessments of bioimage analysis solutions are often arduous, namely due to limited availability of public bioimage D.V. Sorokin was with the Centre for Biomedical Image Analysis, Faculty of Informatics, Masaryk University, 60200 Brno, Czech Republic. He is currently with the Laboratory of Mathematical Methods of Image Processing, Faculty of Computational Mathematics and Cybernetics, Lomonosov Moscow State University, 119991 Moscow, Russia. He is the corresponding author of this manuscript (e-mail: dsorokin@cs.msu.ru).

I. Peterlík is with the Institute of Computer Science, Masaryk University, 60200 Brno, Czech Republic, and with INRIA, France.

V. Ulman was with the Centre for Biomedical Image Analysis, Faculty of Informatics, Masaryk University, 60200 Brno, Czech Republic. He is currently with the Max Planck Institute of Molecular Cell Biology and Genetics (MPI-CBG), 01307 Dresden, Germany. D. Svoboda, T. Nečasová, K. Morgaenko, L. Eiselleová, L. Tesařová, and M. Maška are with the Centre for Biomedical Image Analysis, Faculty of Informatics, Masaryk University, 60200 Brno, Czech Republic. datasets accompanied by reference annotations [START_REF] Kozubek | Challenges and benchmarks in bioimage analysis[END_REF], the preparation of which is known to be subjective, highly error-prone, and extremely laborious, especially when it comes to manually annotating 3D time-lapse image data by human experts [START_REF] Coutu | Probing cellular processes by long-term live imaging-historic problems and current solutions[END_REF].

Therefore, recent benchmarks of time-lapse bioimage analysis algorithms [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF]- [START_REF] Ulman | An objective comparison of cell-tracking algorithms[END_REF] profit from using computer-simulated image datasets, accompanied by inherently generated reference annotations. The simulation platforms adopted for these benchmarks are capable of generating realistic time-lapse sequences of intracellular and intercellular particles [START_REF] Chenouard | Objective comparison of particle tracking methods[END_REF], and motile cells with nuclear staining [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF]. With the primary aim of stimulating the development of fully 3D routines for filopodium segmentation and tracking [START_REF] Tsygankov | CellGeo: A computational platform for the analysis of shape changes in cells with complex geometries[END_REF]- [START_REF] Jacquemet | FiloQuant reveals increased filopodia density during breast cancer progression[END_REF], a prototype of the simulation system capable of generating 3D time-lapse sequences of single motile cells with filopodial protrusions emerged recently [START_REF] Sorokin | Model-based generation of synthetic 3D time-lapse sequences of motile cells with growing filopodia[END_REF]. However, it can only handle a predefined number of curvilinear filopodia with the whole-simulation-period lifetimes, not being capable of mimicking filopodium appearance, disappearance, and branching events, natural for various cell phenotypes [START_REF] Maška | Automatic quantification of filopodia-based cell migration[END_REF].

In this paper, we considerably extend the above-mentioned simulation system prototype [START_REF] Sorokin | Model-based generation of synthetic 3D time-lapse sequences of motile cells with growing filopodia[END_REF] to be capable of generating realistic 3D time-lapse sequences of single motile cells with filopodial protrusions of user-controlled structural and temporal attributes. The proposed extended simulation system allows one to mimic natural scenarios in which cells asynchronously protrude time-varying numbers of filopodia of different thicknesses, lengths, and levels of branching in response to stimuli from their surrounding environments. It consists of three globally synchronized modules, responsible for stochastic simulations of filopodial evolution [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF], simultaneous deformations of the entire cell with filopodia using linear elasticity [START_REF] Felippa | A systematic approach to the element independent corotational dynamics of finite elements[END_REF], and the synthesis of realistic, time-coherent cell texture utilizing clouds of simulated fluorescent particles in the spatio-temporal domain [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF], the content of which is periodically submitted to a virtual microscope [START_REF] Ulman | Virtual cell imaging: A review on simulation methods employed in image cytometry[END_REF]. The proposed simulation system is validated by generating three 3D time-lapse sequences of lung cancer cells of two different phenotypes. We show that the generated sequences qualitatively and quantitatively resemble their real counterparts acquired using a confocal microscope. Note that the output of the proposed simulation system consists of a 3D synthetic time-lapse sequence and the corresponding labeled masks of the cell body and filopodial branches. Such a reference annotation can directly be used when benchmarking filopodium segmentation and tracking methods, and also when studying to what extent segmentation and tracking errors affect Algorithm 1 The pseudocode of the whole simulation process. Inputs: global time step (∆T ), time limit (T lim ), initial parameters of the model (par). Outputs: annotated 3D time-lapse sequence (frame i for i = 0, . . . , T lim /∆T ).

1: [cell 0 , fList 0 , texture 0 ] ← initializeModel(par) 2: frame 0 ← renderImg(cell 0 , fList 0 , texture 0 ) 3: for (i = 0; i • ∆T <= T lim ; i + +) do 4:

fList new ← growFilo(fList i , ∆T ) 5: [cell i+1 , fList i+1 ] ← deformModel(cell i , fList new , ∆T ) 6:
texture i+1 ← updTexture(cell i , cell i+1 , fList i , fList i+1 , texture i , ∆T )

7:
frame i+1 ← renderImg(cell i+1 , fList i+1 , texture i+1 ) 8: end for the extraction of quantitative parameters of crawling cells with filopodial protrusions, such as the number, length, and lifetime of filopodia, or their protrusion and retraction rates.

The rest of this paper is organized as follows. The proposed simulation system is described in Section II. Its experimental validation can be found in Section III. Concluding remarks and suggestions for future work are finally given in Section IV.

II. METHODS OF THE SIMULATION SYSTEM

The simulation of a 3D time-lapse sequence of one motile cell with filopodial protrusions consists of three isolated parts: modeling of filopodium growth and branching (Section II-A), simultaneous elastic deformation of the cell body and individual filopodia (Section II-B), and finally generation of the realistic cell texture (Section II-C). Each part is addressed by an independent software component: geometry module (GM), deformation module (DM), and texture module (TM). The three modules are executed from the main component, a supervisor, that primarily controls the data exchange between the modules, handles their parameters, and synchronizes the time.

The whole simulation process (see Alg. 1) starts by initializing the cell body (cell 0 ) and individual filopodia (fList 0 ), setting their growth parameters, initial lengths l 0 , and the positions where they are attached to the cell body, generating the initial cell texture (texture 0 ), and finally rendering the first frame of the time-lapse sequence (frame 0 ). In our data representation, each cell body, cell i , includes only its geometry stored as a tetrahedral mesh, whereas each filopodium in fList i is represented by a tree structure that contains the data related to the filopodium geometry and evolution (see Section II-A4 for more details).

After the initialization, the simulation process continues iteratively. In each iteration, first, new filopodia may appear or some of the existing ones may disappear, depending on predefined probabilities P a and P d . Filopodium appearance is implemented by adding a new element to fList i and running the initialization procedure for it. Filopodium disappearance is implemented by the 5-fold increase of its retrograde flow (see Section II-A for more details), leading to a natural retraction of the filopodium into the cell body. Second, the morphology of individual filopodia is evolved (fList new ), which allows their branching or the changes in their lengths and widths. Third, the cell body (cell i ) and the evolved filopodia are deformed using pressures applied to the cell body and forces applied to the tips of filopodial branches, yielding the deformed cell body (cell i+1 ) and the deformed filopodia (fList i+1 ). In case a filopodial branch has protruded, a new tree node is appended to its tip. In case it has shortened, the corresponding tree node is removed. After having the geometry of the whole cell updated, the cell texture must be recomputed accordingly, which boils down to updating spatial positions of the simulated fluorescent particles. Finally, the current frame of the simulated sequence (frame i+1 ) is created by rendering the model.

The modules work with their own internal time steps, regularly exporting their particular outputs. A user-defined global time interval ∆T , during which each of the three modules is operating, corresponds to the time step of the generated timelapse sequence. The simulation process is terminated after the total accumulated time reaches a pre-defined limit T lim ≥ ∆T .

The labeled masks of the cell body and filopodial branches, establishing the ground truth for the generated 3D time-lapse sequence, are formed by projecting and rasterizing the meshes of the cell body and filopodial branches to the image domain.

A. Modeling of filopodium growth and branching

In this section, the mathematical models, data structures, and algorithms employed in the geometry module are described.

1) Filopodium evolution on molecular level: A filopodium is a bundle of n f actin filaments enclosed by the cell membrane [START_REF] Carlier | Actin-based Motility[END_REF] (see Fig. 1). Its length is defined as the length of the longest actin filament l = max(l n ) where n is the actin filament index. Each actin filament consists of actin monomers connected with a half-monomer-size overlap of 2.7 nm [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF]. The evolution of actin filaments is affected by four processes: diffusion of free actin monomers along the filopodium, polymerization and depolymerization of actin filaments, fluctuation of the membrane under load, slowing down the polymerization rate of actin filaments, and a constantvelocity retrograde flow that continuously pulls actin filaments into the cell body, modeled on a molecular level using the Gillespie algorithm [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF].

The filopodium volume is discretized by dividing it into J compartments, with the compartment length of l c = 50 nm. Let a j , j = 1, . . . , J, denote the number of actin monomers in the j-th compartment, being indexed from the place where the filopodium attaches to the cell body (filopodium base). The diffusion of each free actin monomer is defined as 1D random walk [START_REF] Erban | Reactive boundary conditions for stochastic simulations of reaction-diffusion processes[END_REF] with reflective boundary conditions at the filopodium tips (P 1→0 = 0). At the filopodium base, the number of actin monomers remains unchanged at the constant bulk value (a 0 = 5) because the number of free actin monomers in the cell body is significantly higher than in the filopodium. Since free actin monomers move between neighboring compartments and stay in the same compartment with equal probability, the reaction time of free actin monomer diffusion is defined as ∆t = l 2 c 3D , with D being the diffusion constant [START_REF] Erban | Reactive boundary conditions for stochastic simulations of reaction-diffusion processes[END_REF].

The actin filament polymerization rate is given [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF] as

k n = k 0 exp (-f n δ/(k B T ))
where k 0 is the "bare" filament polymerization rate and the exponential term accounts for the probability of opening the gap larger than δ between the tip of the n-th filament and the fluctuating membrane. The term k B T represents the thermal energy and f n is the membrane loading force defined for every filament depending on membrane fluctuations [START_REF] Lan | The stochastic dynamics of filopodial growth[END_REF]:

f n = erfc(l -l n , σ) i erfc(l -l i , σ) f
where σ is the average membrane fluctuation amplitude and f is the total membrane force. The actin filament depolymerization rate does not depend on the membrane loading force and has the same value k d for all filaments. The motion equation for each filament length l n is given as:

∆l n = -ν retr ∆t + ξδ -ηδ
where ∆l n is the change of filament length over the reaction time ∆t, ν retr is the retrograde flow velocity, and ξ and η are the binary indicators of polymerization and depolymerization events, respectively. Every iteration of the Gillespie simulation involves three steps [START_REF] Erban | A practical guide to stochastic simulations of reaction-diffusion processes[END_REF]. First, two random numbers, r 1 and r 2 , uniformly distributed in the [0, 1] interval are generated. Next, assuming k d ≤ k n without loss of generality, the type of event to occur in a particular iteration is established as depolymerization

(ξ = 0; η = 1) if 0 ≤ r 1 < k d k d +kn , polymerization (ξ = 1; η = 0) if k d k d +kn ≤ r 1 < kn k d +kn , or free actin monomer diffusion (ξ = η = 0) if r 1 ≥ kn k d +kn .
Finally, depending on the event type, the number of free actin monomers in each compartment is updated, and the reaction time ∆t is calculated as either ∆t = l 2 c /(3D) in case of diffusion or ∆t = (a tip n (k d + k n )) -1 ln(1/r 2 ) in case of polymerization or depolymerization, with a tip n being the number of free actin monomers in the filopodium compartment that corresponds to the tip of the n-th filament, not necessarily being the first one. The Gillespie simulations are performed individually for each actin filament over the global time step ∆T .

2) Filopodium branching: The filaments in the bundle are tied together by special binding proteins. In case the proteins are not up to tie the complete bundle together, the filopodium branches. Without loss of generality, we assume that it happens at the interface of two successive compartments (see Fig. 1). The new branch volume is again discretized by dividing it into the compartments of the length l c , each having its own number a j of free actin monomers. The values of a j are changing over the Gillespie simulations for the filaments that are located in the corresponding compartments of the new branch.

The branching events depend on the pre-defined probability P b . In our simulation system, every new branch contains a random portion of the filaments from the mother branch, and the rest of the filaments remain in the mother branch. The angle between the new branch and the mother branch is given by a random variable with the uniform distribution U(50 • , 70 • ).

3) Filopodium shape: In order to model the 3D shape of a filopodium, every filopodium tree node has a radius assigned, which linearly decreases from the base point to the tips. The radius at the filopodium base is represented by a random variable with uniform distribution U(r m -r ε , r m + r ε ), where r m and r ε are pre-defined parameters. The radii at the filopodial tips correspond to a third of that at the filopodial base.

4) Data structures: Every filopodium is represented by a tree structure composed of J + 1 nodes, one for each compartment and the root that corresponds to the filopodium base. Each node contains the information related to the filopodium geometry (3D coordinate of the filopodium base for the root node, of the filopodial tip for a leaf node, and of the end of the corresponding compartment for the other nodes; radius at that spatial coordinate; branch ID; and force applied to the node) and to the filopodium evolution (the number of filaments and the concentration of free actin monomers), as shown in Fig. 2. Furthermore, for every filopodium, a list of n f elements, one for each actin filament of the filopodium, is maintained. Each element contains the information about the filament length and the index of the compartment that contains the filament tip.

5) Overall simulation of the filopodial evolution: For each filopodium, given by the corresponding tree structure (fTree) and list of filaments (flm), the simulation process (see Alg. 2) runs with a non-uniform, a priori unknown time step due to a possible occurrence of different events (i.e., diffusion, polymerization, or depolymerization) in the Gillespie algorithm. Each iteration starts by checking the occurrence of filopodium branching, being performed at most once per the global time step ∆T . Next, each filament (f n ) of the filopodium is evolved using the Gillespie algorithm and its updated version (f new n ) is propagated back to the corresponding data structures (fTree and flm). The iterative process finishes after the internal time Algorithm 2 The pseudocode of the geometry module. Inputs: global time step (∆T ), initial filopodium tree structure (fTree), list of n f filaments (flm), branching probability (P b ), geometry module parameters (par GM ). Outputs: updated filopodium tree structure (fTree) and list of filaments (flm).

1: t ← 0 2: branched ← false 3: while (t < ∆T ) do 4:

if not branched and shouldBranch(P b ) then 5:

branchFilopodium(fTree, flm)

6:
branched ← true for

(n = 0; n < n f ; n + +) do 9: [f n ] ← getFilamentTrace(fTree, flm, n) 10: [f new n , ∆t n ] ← evolveFilament(f n , par GM ) 11:
updateFilamentTrace(fTree, flm, f new n )

12:

end for 13:

t ← t + max n (∆t n )
14: end while 15: updateRadii(fTree) 16: setForces(fTree) (t) reaches the global time step ∆T . Finally, the information about the filopodium radius in each compartment and about the forces applied to the filopodial tips is updated accordingly, thus providing the deformation module with the actual information needed to properly deform the evolved filopodium.

B. Deformation of cell body and filopodia

The motion and deformation of the cell is modeled by the deformation module, which is started just after geometry module every global time step. To mimic realistic deformations, the cell body and the filopodia are modeled as elastic objects coupled together via mechanical constraints. Each filopodium and the cell body has its proper dynamics based on the elastic formulation discretized by the finite element method (FEM). The mechanical coupling between a filopodium base and the cell body propagates the excitations in both directions: forces applied to filopodia influence the motion of the cell body and random pressures applied to the cell body surface result in the motion of filopodia. We first characterize the dynamic formulation of filopodial and cell body motions and then, we describe the mechanical coupling.

1) Mechanical model of filopodia: Each filopodium is modeled as a deformable tree of piecewise linear segments of nonzero radii, represented by linear beam elements based on the Timoshenko formulation [START_REF] Przemieniecki | Matrix structural analysis of substructures[END_REF]. While being computationally efficient, the beam elements are typically used to model threadlike structures which undergo bending and display important force responses when subjected to longitudinal deformations.

Let a filopodium ϕ, represented by a fTree structure, be given as a set of nodes {q 0 , . . . q n }, each being an oriented point in 3D space with three positional and three rotational degrees of freedom (DoF). Its dynamics is governed by the Newton's law of motion relating the body acceleration and forces:

M ϕ q = f ϕ -g ϕ (q, q) + H ϕ λ (1)
where M ϕ is an inertia matrix, q is a vector of filopodium node coordinates, f ϕ is a vector of external forces, and g ϕ is a vector of internal elastic forces. The term H ϕ λ represents the coupling between the filopodium and cell body, explained in more detail in Section II-B3.

The internal forces are represented as a general non-linear function that is linearized in each internal time step of the deformation module:

g ϕ (q + dq, q + d q) ≈ g ϕ (q, q) + ∂g ϕ ∂q dq + ∂g ϕ ∂ q d q. (2)
The derivatives K ϕ = ∂gϕ ∂q and B ϕ = ∂gϕ ∂ q are known as stiffness and damping matrices, respectively. In the actual scenario, the inertia matrix M ϕ is constant and diagonal. Since we employ the co-rotational formulation of each beam element [START_REF] Przemieniecki | Matrix structural analysis of substructures[END_REF], [START_REF] Cotin | New approaches to catheter navigation for interventional radiology simulation[END_REF], the stiffness matrix K ϕ is updated in each internal time step because the rotational components depend on the actual deformation. The stiffness matrix K ϕ determines the physical behavior of the filopodium and is parametrized by Young's modulus E ϕ and Poisson ratio ν ϕ . The damping matrix B ϕ is approximated using Rayleigh stiffness r K and Rayleigh mass r M , as B ϕ = r K K ϕ + r M M ϕ . The vector of external forces f ϕ represents the excitation forces applied to the filopodial tips. For every filopodial branch, the force is applied to its tip in random direction orthogonal to its last segment, and the force magnitude is given as the simulation parameter F tip .

The integration of ( 1) is carried out using the implicit Euler scheme [START_REF] Baraff | Large steps in cloth simulation[END_REF] computed with an internal time step h. Denoting t the previous DM time step, the linear system in the actual DM time step t + h is given [START_REF] Peterlik | Constraint-based haptic rendering of multirate compliant mechanisms[END_REF] as:

(M ϕ -hB ϕ -h 2 K ϕ ) Aϕ d q = hg ϕ -h 2 K ϕ qt +f ϕ bϕ +hH ϕ λ (3)
where M ϕ is constant, and B ϕ , K ϕ , and g ϕ are computed using the position q t obtained in the previous internal time step.

The linear system is solved by the square-root-free Cholesky decomposition in the form of A ϕ = LDL , with L being a unit lower triangular matrix and D being a diagonal matrix.

2) Mechanical model of the cell body: We use the same approach based on the implicit Euler method to model the cell body, however, unlike in the case of filopodia, the cell body is modeled using 3D co-rotational elasticity discretized by linear tetrahedra. Therefore, the cell body domain is given by a set of nodes {p 0 , . . . , p m } where each p i is given by three positional degrees of freedom. The discretized Newton's equation of cell body motion is

M c p = f c -g c (p, ṗ) + H c λ (4)
where M c is the constant inertia matrix, p is a vector of cell body node coordinates, and f c and g c are the external and internal forces of the cell body, respectively. The term H c λ represents the coupling between the cell body and the filopodia, explained in more detail in Section II-B3. The same linearization procedure is applied as in the case of filopodium: given the constant matrix M c and K c , B c , f c and g c computed using the positions p t obtained in the previous DM time step, the dynamic system in the DM actual step t + h is given as

(M c -hB c -h 2 K c ) Ac d ṗ = hg c -h 2 K c ṗt +f c bc +hH c λ. (5)
The global stiffness matrix K c is assembled from local matrices K cτ ; each K cτ computed for a tetrahedron τ is a 4×4 block matrix where each block is of size 3×3 and is determined by Young's modulus E c and Poisson ratio ν c of the cell. Please see [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF] for further details about the co-rotated formulation of linear elasticity and the definition of matrix K c . The external forces f c applied to the cell body result from the random pressures acting on the cell surface [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF]. The magnitude of the random pressures is set as the simulation parameter P cb .

3) Coupled dynamic simulation of cell with filopodia: Without loss of generality, we describe a system composed of a cell body and one filopodium ϕ and it is supposed that initially the 6DoF base q 0 0 of the filopodium is located inside a tetrahedron τ ϕ . During the simulation, the initial relative position and orientation of q 0 0 with respect to τ ϕ are maintained: therefore, we introduce mapped base r ϕ : initially, r ϕ = q 0 0 and during the simulation, its position is updated together according to the actual position of the tetrahedron τ ϕ .

The relative position and the orientation of r ϕ are given by the barycentric coordinates β ϕ of r ϕ with respect to τ ϕ and the rotation matrix Ω ϕ , which are initially precomputed and remain constant during the simulation [START_REF] Peterlík | Modeling and real-time simulation of a vascularized liver tissue[END_REF].

In the actual DM time step t + h, the coupled dynamics of the cell with filopodium is computed as follows:

• Free motion: the vector of constraint forces is set to zero (λ = 0), and (3) and ( 5) are solved for d q and d ṗ, respectively. These are then integrated to obtain positions of both the filopodium and the cell body according to unconstrained dynamics: q free = q t + h( qt + d q), p free = p t + h( ṗt + d ṗ). Since the tetrahedron τ ϕ has moved, it is necessary to update the mapped filopodium base r ϕ using the precomputed barycentric coordinates β ϕ and its orientation using the rotation matrix Ω ϕ . • Computation of coupling forces: Due to the independent dynamics computed by the free motion, the actual filopodium base q 0 and mapped filopodium root r ϕ do not necessarily coincide. However, to maintain the coupling between the two objects, it is necessary to impose q 0 = r ϕ , i. e. δ = |r ϕ -q 0 | = 0. This is done by computing the constraint force λ that must be applied to both objects in order to achieve the equality δ = 0. As shown in [START_REF] Peterlik | Constraint-based haptic rendering of multirate compliant mechanisms[END_REF]

δ = δ free + h[H ϕ A -1 ϕ H ϕ + H c A -1 c H c ]λ (6)
where δ free is a constraint violation after the free motion, and A -1 ϕ and A -1 c are the inverted system matrices computed during the free motion in (3) and [START_REF] Jacquemet | Filopodia in cell adhesion, 3D migration and cancer cell invasion[END_REF]. The matrices H ϕ and H c map the force vector λ to the degrees of freedom of respective objects. In the case of filopodium, H ϕ is a 6 × 6N matrix (N is the number of nodes in the filopodium), which is zero everywhere except for the leftmost block, corresponding to the first node, which is given as the 6×6 identity matrix. The 6 × 3M matrix (M is the number of nodes in the cell body) H c maps the 6DoF forces λ to four linear forces applied to the nodes of the tetrahedron τ ϕ . Therefore, H c contains four non-zero blocks of size 6 × 3: one for each vertex I of the tetrahedron τ ϕ . Each block is given as

H cI =   β I 0 0 0 -∆ z I ∆ y I 0 β I 0 ∆ z I 0 -∆ x I 0 0 β I -∆ y I ∆ x I 0   (7)
where β I are the barycentric coordinates of r ϕ and ∆ I = p I -r ϕ is the line given by the updated mapped position of the filopodium base and the vertex I of the tetrahedron τ ϕ . The system ( 6) is solved using the iterative Gauss-Seidel method. • Corrective motion: The corrected positions of both the cell body and filopodium are obtained by applying the coupling forces λ to the configurations obtained by the free motion. This is done by computing

q t+h = q free + A -1 ϕ H ϕ λ (8)
p t+h = p free + A -1 c H c λ. (9)
4) Cell body motion: At the beginning of the execution of the deformation module on every global time step, the cell centroid coordinates are moved to a new location and rotated by a random angle in 3D to mimic the random motion of the cell. The cell motion is incorporated to the model by applying the force to the center of the cell body (for translation), and on several nodes close to the center of the cell body (for rotation). Both translational and rotational forces are included in the external forces f c in (5) thus the other parts of the cell are influenced taking into account the elasticity of the cell which results in motion that resembles amoeboid-like motility. These forces are defined using a trajectory generated beforehand which is a simulation parameter. The cell rotation is defined by the random rotation axis and the random rotation angle γ, which is also a simulation parameter.

C. Generation of cell texture

The cell texture is treated with three main objectives in mind: high realism of the texture, ability to adapt to cell shape changes, and visual coherence between consecutive frames of the generated time-lapse sequence. Two different paradigms are employed. For the cell body, which is not changing its volume rapidly, a texture visual representation is properly initiated for the first time point, and is updated according to tracked changes in cell geometries in consecutive time points. For every filopodial branch, which besides bending is also growing or retracting, a new texture visual representation is created for every time point.

The visual representation of the simulated cell is implemented using a set of numerous texture particles where each particle represents a clump of fluorescent molecules [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF], [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF], [START_REF] Ghaye | Simulated biological cells for receptor counting in fluorescence imaging[END_REF]. This is a consequence of the fluorescence microscopy principle that only cell parts stained with a particular fluorescent dye (GFP, Hoechst, etc.) are visible in acquired images. In this study, we mimic transient cell transfection with the green fluorescent protein conjugated to actin. Such a labeling, required for a good definition of the filopodial structures, produces a heterogeneous fluorescent signal within the cell body volume, being concentrated in the cortical cytoskeleton of the cell and leaving the inner parts of the cell practically unstained.

Every texture particle is represented by a floating-point micrometer-scaled three-dimensional coordinate vector, and initiated with a fixed amount of fluorescence intensity it carries. Projecting all particles together into a 3D lattice (of generally anisotropic resolution) and accumulating their fluorescence gives rise to a digital phantom image I : Z 3 → N 0 (see Fig. 3).

To address different geometry developments and labeling outcomes for cell body and for filopodial branches, we treat the ensemble of texture particles in each case separately.

1) Cell body: To initiate the set of texture particles, a cell body mask is created, filled with a Perlin noise [START_REF] Perlin | An image synthesizer[END_REF], and particles are sampled from this. To create the mask, the volumetric tetrahedral mesh from DM is superimposed over the 3D lattice to give rise to a cell body mask image M b : Z 3 → {0, 1} which 1 is set to voxels that intersect with the volume of the mesh, and 0 otherwise.

Perlin noise [START_REF] Perlin | An image synthesizer[END_REF] yields a coherent procedural texture image ρ : Z 3 → R controlled by four parameters influencing its coarseness C, smoothness S, flickering F l, and fineness F i. The coarseness C is a number controlling the magnification of the blobs in the texture and we introduce it to make our model scalable. The remaining three parameters are introduced in [START_REF] Perlin | An image synthesizer[END_REF]. The Perlin texture is composed of periodic functions, i.e., the domain of texture is unlimited. In order to always get a unique texture image we have to apply a random coordinate shift S[.] on the generated texture ρ. This type of texture images is commonly accepted as a good approximation of randomized cellular structures [START_REF] Rousselle | Chromatin texture analysis in living cells[END_REF].

A mosaic of Perlin noises is used to model the initial texture of a cell body. Inner part of the cell body has very fine and dark structure. This part is generated using a fine-grained Perlin noise ρ ci (parameters C ci , S ci , F l ci , and F i ci ):

Ĩci = S[ρ ci ]. ( 10
)
Compared to the inner structure, the boundary of the cell body (which corresponds to the cortical part) is much brighter and coarser. Here, we employ a coarse-grained Perlin noise ρ cb (parameters C cb , S cb , F l cb , and F i cb ):

Ĩcb = S[ρ cb ]. (11)
To create a texture image Ĩ to initialize texture particles with, a linear blending between Ĩcb and Ĩci is performed and the product is masked with M b . Let w be the width of the cell body's perimeter volume that represents the cell boundary with coarser texture, let s be the width of a transition strip between the two textures, and EDM be the Euclidean distance map of the foreground voxels in M b . The texture image is then defined as:

Ĩ(x) = α(x) • Ĩci (x) + (1 -α(x)) • Ĩcb (x) • M b (x) , (12)
where α(x) = min (s, max(0, EDM(x) -w)) /s (in this work we used w = 0.1 µm and s = 0.05 µm). To obtain an initial distribution of particles in the cell body, the particles are introduced for every texture voxel of the image Ĩ in a number that is proportional to that voxel intensity. The set of all cell body particles' coordinates is marked as {p i } and its size is |{p i }| = x Ĩ(x)/Q , with Q being a pre-defined quantization factor. The number |{p i }| is constant throughout the entire simulation.

The geometry of the cell body is represented with volumetric tetrahedral mesh. The GM provides initial body mesh which is further deformed by DM every global time step based on the elastic behavior of the cell. This way, the DM inherently generates the displacements in all mesh vertices. By interpolating the displacements between the mesh vertices one can form a vector field that represents geometry change from the previous to the current time point. The coordinates pi of the texture particles, fitting originally to the geometry for the previous time point, are updated according to this vector field. We have previously shown that applying this principle to the texture representation satisfies the three objectives [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF] above. Since the deformation field generated by the DM is smooth due to the chosen elasticity model, the update of texture cannot introduce abrupt changes in voxel intensities. However, the condensation or spread of the texture, leading to a local increase or decrease in voxel intensities, can occur in case the corresponding part of the cell volume changes, for example nearby protruding cell membrane.

2) Filopodial branch: Geometry of each filopodial branch ϕ is fed into the TM in the form of a Timoshenko beams, that is, as a list of points q i that sample the beams, with radius r i associated to every point. This geometry can also be understood as a piece-wise conical object that poses a base-totip central axis and whose width decreases from the base to the tip. Indeed, a surface mesh representing the filopodial branch is built and rasterized (see Fig. 3) into a mask image M ϕ : Z 3 → {0, 1} in which 1 is set to voxels that are inside this mesh, and 0 otherwise. Considering any conical segment of the filopodium, e.g., between points q k and q k+1 , two circles are constructed, sampled with 10 points and surface triangles are created between two consecutive points on one circle and one point on the other circle. A circle centered at point q k is of radius r k and coincides with a plane with a normal vector q k+1 -q k-1 . The normal vectors used for circles associated with the first and last points of the filopodium branch are colinear with the first and last segments, respectively.

The coordinates qi of texture particles representing a filopodium ϕ are introduced in the similar fashion. In real image data, the filopodia are visually recognized as bright compartments without an apparent internal structure. To achieve the same visual experience, the volume of a filopodium is populated with a larger number of uniformly distributed particles, making the resulting texture insensitive to particular precise coordinates of the individual particles. Any beam is, therefore, regularly sampled in |q k -q k+1 |/M axial steps (see Table I for the parameter M axial value), and a disc centered at every sampling point is assumed. The disc is further assumed perpendicular to the axis at the point, and with radius r linearly interpolated between r k and r k+1 . Finally, the disc is populated with M radial • r 2 particles, where M radial is a constant influencing the density of particles in this disc (see Table I). Iterating over all sampling points and over all beams, the set of all filopodial branch particles' coordinates {q i } is created. Due to the changes in the geometry and volume of the filopodium, the number of created particles usually changes between time steps, but their density is maintained. This allows the texture particles to be always created from scratch to the current filopodium geometry in every time point without compromising visual coherency of the filopodium texture in successive time points.

For any time point, it now remains to construct output images. The labeled mask is obtained as a weighted sum of all the mutually disjoint masks M b and M ϕ , where the weights are the expected numerical labels of individual components. The phantom image I(x) is constructed using union {p i } ∪ {q i } of the sets with texture particles' coordinates that existed at the respective time point:

I(x) = ζ∈{pi}∪{qi} χ(ζ, x) (13)
where

χ(ζ, x) = Q if particle ζ falls into voxel x, otherwise χ(ζ, x) = 0.

D. Simulation of optical system and acquisition device

To obtain the final image that is visually similar to a real microscopy one, the phantom image I is submitted to a virtual microscope [START_REF] Svoboda | Generation of digital phantoms of cell nuclei and simulation of image formation in 3D image cytometry[END_REF]. In particular, the phantom image is first blurred by an experimental point spread function, measured on a real optical system, and second affected by an uneven illumination phenomenon [START_REF] Svoboda | On proper simulation of phenomena influencing image formation in fluorescence microscopy[END_REF]. Next, it is submitted to a virtual CCD camera. The uncertainty of the photo detection in the camera is modeled with photon shot noise that is subject to Poisson distribution. Finally, the image is affected by amplifier read out, and dark current noise, all of them according to the camera documentation (see Fig. 3). As the final step, the anisotropy of the final synthetic image is obtained by subsampling the synthetic image along z-axis.

III. RESULTS

To demonstrate the capability and flexibility of the proposed simulation system, we imitated actin-stained A549 lung cancer cells stably transfected with two genetically different variants of the CRMP-2 protein: overexpressing and phospho-defective, which are the representatives of mesenchymal migration typical for cancer extravasation [START_REF] Jacquemet | Filopodia in cell adhesion, 3D migration and cancer cell invasion[END_REF]. Depending on its phosphorylation status, CRMP-2 is a cytosolic phosphoprotein that partly contributes to the microtubule dynamics by either promoting or impairing the ability of adding tubulin heterodimers to the microtubule plus ends. Therefore, it indirectly affects filopodial formation via microtubule scaffolding that bidirectionally interacts with actin filaments, both actively regulating the cell cytoskeleton shape [START_REF] Mattila | Filopodia: Molecular architecture and cellular functions[END_REF].

As a reference, we considered real image data with the voxel size of 0.126×0.126×1.0 µm, being acquired every two minutes over one hour using a spinning disk confocal microscope equipped with a water Plan-Apochromatic 63×/1.20 objective lens [START_REF] Maška | Automatic quantification of filopodia-based cell migration[END_REF]. Whereas the overexpressing cells had variable cell body texture and protruded very short, single-branch filopodia of low motility that frequently appeared and disappeared, the phospho-defective cells had less heterogeneous cell body texture and protruded long, branching filopodia of high motility with whole-experiment lifetimes. As a counterpart to three real 3D time-lapse sequences of each of the two cell phenotypes, we generated three synthetic 3D time-lapse sequences of either cell phenotype using the proposed simulation system. In total, we analyzed 180 real and 180 synthetic images to qualitatively and quantitatively evaluate the similarity between the real and computer-generated image data.

C ci U (1, 2) 1 S ci 0.3 0.3 F l ci 1 1 F i ci 6 6 C cb 2.5 2.5 S cb 0.2 0.2 F l cb 1 1 F i cb 6 6 Q 64 64 M axial 10 5 M radial 12 4
The computer-generated image data, including its complete segmentation and tracking annotations, is made publicly available at http://cbia.fi.muni.cz/research/simulations/filogen.html.

A. Simulation parameters

The initial cell bodies were generated using the digital cell phantom simulator [START_REF] Svoboda | On proper simulation of chromatin structure in static images as well as in time-lapse sequences in fluorescence microscopy[END_REF]. The parameter configurations used for generating synthetic 3D time-lapse sequences of both overexpressing and phospho-defective cells are listed in Table I. They reflect the diversity in the initial number (N ϕ ), level of branching (P b ), lifetime (P a and P d ), shape, and dynamics of filopodia, as well as in the heterogeneity of cell body textures (C ci ). In particular, the longer (l 0 ) and thicker (r m and r ε ) filopodia of the phospho-defective cells require less force to be bended (F tip ) than those of the overexpressing cells. The number of actin filaments (n f ) in the filopodia was chosen to be between 10 and 30, not to limit the filopodial growth by the plasma membrane resistance or insufficient free actin diffusion at filopodial tips [START_REF] Mogilner | The physics of filopodial protrusion[END_REF]. The motion of both cell phenotypes was given by mimicking random-walk trajectories of the real cell centroids extracted from the available segmentations of the reference real image data [START_REF] Maška | Automatic quantification of filopodia-based cell migration[END_REF]. In order to imitate a diversity of the cell texture among the individual real cells, the value of C ci parameter was generated for each overexpressing cell individually as a uniformly distributed random number. Two synthetic 3D time-lapse sequences that were generated using the proposed simulation system with the listed parameter configurations are shown in Figs. 4 and5, demonstrating the growth and motion of filopodia of a phospho-defective and an overexpressing cell, respectively. A qualitative comparison of the real and synthetic image data showing overexpressing and phoshodefective cells is given in Fig. 6. It can be seen that the textures of cell bodies and filopodia in the real and synthetic image data are visually comparable.

B. Validation of synthetic image data

The generated synthetic time-lapse datasets were thoroughly designed to aptly mimic the shape, size, internal structure, and filopodial tip dynamics of the reference real cells. To validate that our simulation system produces plausible image data, we performed several tests that quantitatively assess the similarity between the synthetic datasets and their real counterparts.

1) Texture of cell body: In order to avoid the influence of the background noise, we calculated gray level co-occurrence matrices over the cell body interior regions (31×31×3 voxels, see Fig. 6b) and followed the protocol introduced in [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF]. In particular, five Haralick's texture descriptors: Texture contrast, Texture correlation, Texture homogeneity, Maximal correlation coefficient, and Entropy, were computed for all synthetic and real images, and their distributions, at the level of cell phenotypes, were mutually compared. To this end, the Kolmogorov-Smirnov test with quantile-quantile plots were used for a nonparametric comparison of the underlying distributions obtained for the real and synthetic image data (see Fig. 7).

The visual inspection of the quantile-quantile plots depicted in Fig. 7 reveals a high level of similarity because the points fit the straight line, representing the ideal identical distributions, most of the range. Except for one single case, all Kolmogorov-Smirnov tests retained the null hypothesis (p > 0.05), stating that the distribution of texture descriptors is the same in both groups. Note that the level of similarity reported here is higher than that reported in the previous study [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF].

Despite an apparent texture heterogeneity in the real overexpressing cells, the results confirmed the ability of our system to imitate real cell body textures with comparable distributions of the evaluated texture descriptors. These descriptors computed for the synthetic phospho-defective cells were also comparable to the real ones, only their Entropy was slightly shifted to the lower values when compared to the distribution extracted from the real image data.

2) Motion of filopodial tips: To obtain the reference motion patterns of filopodial tips for each mimicked cell phenotype, three experts manually annotated spatio-temporal positions of individual filopodial tips available in the real image data. Next, the three created manual annotations were combined together by maximizing the overall Linear Oriented Forests Matching measure [START_REF] Maška | Particle tracking accuracy measurement based on comparison of linear oriented forests[END_REF] between the combined, reference annotation and each of the three manual annotations. Finally, each reference trajectory (i.e., the trajectory of one particular filopodial tip in the reference annotation) was analyzed as described in [START_REF] Wu | Statistical analysis of cell migration in 3D using the anisotropic persistent random walk model[END_REF]. In particular, we studied the orientations of the filopodial tip movement between two successive frames with respect to the primary trajectory direction. The possible orientations, always qq q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q qq q q qq q q q q q q q q q q q q q q q q q q q q q q q 2.97 Texture contrast p = 0.400 qq q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q 0.00 0.12 0.23 0.35 0.00 0.12 0.23

0.35

Real image data Synthetic image data Texture correlation p = 0.164 q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q 0.32 Texture homogeneity p = 0.948 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q qqq q -0.01 0.04 0.09 0.14 -0.01 0.04 0.09 0.14

Real image data Synthetic image data Max. correl. coefficient p = 0.515 q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q 2.08 q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q 3.26 Texture contrast p = 0.948 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q 0.00 0.06 0.12 0.18 0.00 0.06 0.12 0.18

Real image data Synthetic image data Texture correlation p = 0.081 q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.32 Texture homogeneity p = 0.081 q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q -0. Real image data Synthetic image data Max. correl. coefficient p = 0.872 q q qq qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2.07 falling in the [0, 180) interval, were binned into 12 subintervals with the width of 15 degrees (see Fig. 8). The reference motion patterns of filopodial tips for the overexpressing and phosphodefective cells were established by aggregating total distances within individual angular subintervals across all filopodial tips of either cell phenotype and by normalizing them to produce 1 when summed over all angular subintervals. To assess the similarity of motion patterns of filopodial tips displayed in the synthetic image data, we applied the identical routine on its inherently generated reference annotations, and calculated the histogram intersection measure [START_REF] Cha | On measuring the distance between histograms[END_REF] between the motion patterns extracted from the real and synthetic image data of the particular cell phenotype (see Fig. 8). By inspecting the content of Fig. 8, it can be concluded that the proposed Fig. 8. The distance percentage traveled by all filopodial tips of overexpressing (upper) and phospho-defective cells (lower) in the real and synthetic image data at different orientations with respect to the primary trajectory direction of individual filopodial tips. The histograms are accompanied by the information on their similarity calculated using the histogram intersection measure.

simulation system is capable of plausibly mimicking random motion patterns of real filopodial tips, with the high degree of similarity exceeding the level of 0.87 on the [0, 1] scale.

C. Implementation details and execution time

The three modules and the supervisor involved our simsystem were created as separate software components that exchange their data by and writing to files. GM and the supervisor implemented in Matlab R2017a. DM and TM were implemented in C++, the former being built on the SOFA framework [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF]. The modules and the supervisor were integrated into a pre-configured virtual machine, which is made publicly available, together with the source codes, at http://cbia.fi.muni.cz/research/simulations/filogen.html.

To provide the reader with an estimate of how long the generation of the synthetic image data took, we measured the execution time of each phase involved in the simulation process on a common workstation (Intel Xeon QuadCore 2.83 GHz, 32 GB RAM, Ubuntu 16.04 LTS).

The measured average execution times for generating one synthetic 3D image of a single overexpressing and phosphodefective cell are listed in Table II. By comparison with 346 seconds spent on average by one human expert when manually annotating one randomly selected slice for segmentation and the whole 3D image for filopodial tip tracking, it can be concluded that the proposed simulation system substantially reduces the time needed to establish a benchmarking dataset of completely annotated image data for objective comparison of filopodium segmentation and tracking methods.

IV. CONCLUDING REMARKS

We have presented a novel approach to generating synthetic 3D time-lapse sequences of single motile cells with filopodial protrusions of flexible structural and temporal attributes. Being composed of three globally synchronized modules, our system is capable of simulating the filopodium growth and morphology changes, modeling the joint motion and deformation of the cell body with a time-varying number of filopodia, and generating the realistic texture of the entire cell together with mimicking the process of fluorescence microscopy image formation. The filopodium evolution is modeled using stochastic simulations of filamentous growth on a molecular level. The joint motion and deformation of the cell body with multiple filopodia is built on the theory of elasticity and is efficiently implemented using the FEM. The realistic time-coherent cell texture is synthesized using the simulated fluorescent particles in the spatio-temporal domain, being periodically submitted to a virtual microscope. The proposed simulation system was validated by comparing three generated 3D time-lapse sequences of single overexpressing and phospho-defective cells to their real counterparts. The highly flexible modules involved in our system allowed us to fine-tune the synthetic image data in a way to achieve high structural, visual, and quantitative similarity to the real image data in terms of cell body texture and filopodial tip motion.

Such a level of similarity to the real image data and the fact that our simulation system inherently generates completely annotated image data have stimulated the development of fully 3D bioimage analysis pipelines for filopodium segmentation and tracking [START_REF] Castilla | Three-dimensional lifetime quantification of filopodia in migrating cancer cells[END_REF], including those based on nowadays extremely popular machine-learning approaches that often suffer from a limited amount of annotated training data. Furthermore, the generation of such a synthetic dataset is significantly faster than the manual annotation of the real image data. Last but not least, our system is capable of generating synthetic image data of different spatio-temporal resolutions, which is something very difficult to achieve in real microscopy setups due to their physics-based limitations. Such a feature may help in objective benchmarking of the developed bioimage analysis algorithms and in deeper understanding of their limitations, thus providing users as well as algorithm developers with practical clues on the algorithm behavior.

In spite of being flexible and mathematically sound, our system has its limitations that stem from the design of individual modules. This is particularly true for GM and DM. The former has been designed to simulate the evolution of filopodia, barely allowing straightforward modifications that would facilitate to model the evolution of other types of cell protrusions, such as blebs. Furthermore, the current model of filopodial branching, being built on a random process, could possibly be improved by modeling the dynamics of cross-linking proteins to control the bundling and deviation of actin filaments in a biologically more relevant fashion [START_REF] Aratyn | Intrinsic dynamic behavior of fascin in filopodia[END_REF]. The latter module strictly requires the modeled cell geometries to be decomposable into a reasonable number of volumetric elements that are supported by the SOFA framework [START_REF] Faure | Sofa: A multi-model framework for interactive physical simulation[END_REF], and to behave according to the theory of elasticity that forms widely acceptable principles of plasma membrane dynamics of diverse cell phenotypes [START_REF] Danuser | Mathematical modeling of eukaryotic cell migration: Insights beyond experiments[END_REF]. Thus, the requirement for such elastic behavior may allow one to mimic sophisticated cell migration modes, such as mesenchymal and amoeboid motility [START_REF] Driscoll | Quantifying modes of 3D cell migration[END_REF]. Finally, concerning TM, the clouds of simulated fluorescent particles can represent a broad range of fluorescent reporters [START_REF] Svoboda | Mitogen: A framework for generating 3D synthetic time-lapse sequences of cell populations in fluorescence microscopy[END_REF]. Furthermore, their spatial resolution is not bound to any particular image grid, which allows one to interconnect this texture model with spatially high-resolution models of cell shape development [START_REF] Shao | Coupling actin flow, adhesion, and morphology in a computational cell motility model[END_REF]- [START_REF] Prost | Active gel physics[END_REF].

In our future work, we intend to further extend and adapt the powerful mathematical models employed in all three modules of our simulation system, with the primary aim of simulating scenarios with multiple, mutually interacting cells with filopodial protrusions, which exhibit different cell migration modes and control filopodial branching via insufficient concentrations of cross-linking proteins within the bundles of actin filaments.

Fig. 1 .

 1 Fig. 1. An illustration of the molecular model of an evolving filopodium, with the compartment size lc and the filament length ln. Note that the numbers of filaments and of free actin monomers were purposely decreased compared to the real model values to make the figure content visually simple and intuitive.

Fig. 2 .

 2 Fig. 2. The fTree structure representing the purely illustrative model of the filopodium in Fig. 1. The upper part of each node contains the information for geometry modeling and deformation, the lower part contains the information for filopodium evolution.

Fig. 3 .

 3 Fig. 3. Illustration of the data processed during the pipeline of the TM module: Surface of an input mesh, rasterized as a labeled mask, obtained phantom image and final output image. All images are z-axis maximum projections and are enhanced for visual purposes. The voxel size is 0.126 × 0.126 µm.

Fig. 5 .

 5 Fig. 5. Two maximum intensity projections of a synthetic 3D time-lapse sequence of an overexpressing cell. The voxel size is 0.126 × 0.126 × 1.0 µm.

Fig. 4 .

 4 Fig. 4. Two maximum intensity projections of a synthetic 3D time-lapse sequence of a phospho-defective cell. The voxel size is 0.126 × 0.126 × 1.0 µm.

Fig. 6 .

 6 Fig. 6. A 3D visualization of real (a, c) and synthetic cells (b, d) using XY and XZ slices. The overexpressing cells (a, b) have rather short filopodia, whereas the phospho-defective cells (c, d) have long filopodia with several branches. The overexpressing cells are visualized using single slices (green ticks show the selected slices). The phospho-defective cells are visualized using maximum intensity projections as the filopodia are long and difficult to see in a single slice. The green box in Fig. 6b depicts a sample 31×31×3 region of interest used for texture analysis.

Fig. 7 .

 7 Fig.7. Quantile-quantile plots and boxplots for the distributions of Haralick's textural descriptors computed over the cell body interior regions (see Fig.6b) cropped from the real (R) and synthetic (S) image data of overexpressing (left column) and phospho-defective (right column) cells. The displayed quantiles are accompanied by the ideal identical distribution (straight line) and p-values of the Kolmogorov-Smirnov test. The whiskers of boxplots correspond to the minimum and maximum values, the boxes represent the interquartile ranges, and the thick marks stand for the median values.

TABLE I PARAMETER

 I VALUES USED FOR GENERATING SYNTHETIC 3D TIME-LAPSE SEQUENCES OF OVEREXPRESSING AND PHOSPHO-DEFECTIVE CELLS.

	Group	Name Overexpres. Phospho-def. Unit
	Global	∆T	20	20	s
		T lim	580	580	s
		Nϕ	12	5
	Shape	P b	0	0.5
		Pa	0.5	0
		P d	0.5	0
		rm	0.25	1	µm
		rε	0.05	0.3	µm
	Gillespie	D	5	5	µm 2 s -1
		n f	16	16
		l 0	0.081	U (3, 4) µm
		k d k 0	1.4 21.8	1 52.4	s -1 s -1
		k B T	4.1	4.1	pN • nm
		σ	10	10	nm
		f	10	U (1.5, 2.5) pN
		νretr	U (7, 13) U (30, 50) nm/s
	Deformation F tip	25	12.5	µN
	and motion P cb	25	10	kPa
		γ	0.1	0.1	rad
		Eϕ	1	1	kPa
		Ec	50	50	kPa
		νϕ	0.45	0.45
		νc	0.4	0.4
	Texture

TABLE II THE

 II AVERAGE EXECUTION TIME (IN SECONDS) OF EACH PHASE INVOLVED IN THE PROPOSED SIMULATION PROCESS FOR GENERATING ONE 3DIMAGE OF A SINGLE OVEREXPRESSING AND PHOSPHO-DEFECTIVE CELL.

	Phase	Overpressing Phospho-defective
	Geometry module	20	173
	Deformation module	8	3
	Texture module	4	28
	Virtual microscopy	12	36
	Total	44	240

ACKNOWLEDGMENT

The authors would like to thank X. Morales and A. Rouzaut (Center for Applied Medical Research, Spain) for providing the real image data of overexpressing and phospho-defective cells and for their valuable comments on the biological background of the CRMP-2 protein and its mutations.


David Navarre
email: navarre@irit.fr

Philippe Palanque
email: palanque@irit.fr

Arnaud Hamon

Sabrina Della Pasqua

Similarity as a Design Driver for User Interfaces of Dependable Critical Systems

Keywords: UI properties, similarity, dependability, usability, learnability

Assuring that operators will be able to perform their activities even though the interactive system exhibits failures is one of the main issues to address when designing and implementing interactive systems in safety critical contexts. The zero-defect approaches (usually based on formal methods) aim at guaranteeing that the interactive system will be defect free. While this has been proven a good mean for detecting and removing faults and bugs at development time, natural faults (such as bit-flips due to radiations) are beyond their reach. One of the way to tackle this kind of issue is to propose redundant user interfaces offering multiple ways for the user to perform operations. When one of the interaction mean is failing, the operator can select another functional one. However, to avoid errors and increase learnability, it is important to ensure that the various user interfaces are "similar" at presentation and interaction levels. This paper investigates this relation between dependability and similarity for fault-tolerant interactive systems.

Introduction

Usability [START_REF]Ergonomic requirements for office work with visual display terminals (VDT) -Part 11 Guidance on Usability[END_REF] and user experience [START_REF] Hassenzahl | Hedonic and ergonomic quality aspects determine a software's appeal[END_REF] properties have received (and are still receiving) a lot of attention in the area of Human-Computer Interaction to the extent that they are perceived as the main properties to study and consider while designing interactive systems or while performing research activities in HCI. Beyond this main stream of research and design, other more marginal approaches have tried to investigate the relationship between these properties and other ones such as security [START_REF] Sasse | Designing and evaluating usable security and privacy technology[END_REF], accessibility [START_REF]The Road to Accessibility[END_REF][START_REF] Petrie | The relationship between accessibility and usability of websites[END_REF], dependability [START_REF] Fayollas | An Approach for Assessing the Impact of Dependability on Usability: Application to Interactive Cockpits[END_REF] or privacy [START_REF] Gerber | Usability versus privacy instead of usable privacy: Google's balancing act between usability and privacy[END_REF] (among many others).

Each of these specific domains bring specific issues in order to ensure that the associated properties have been taken into account. Taking into account these properties usually requires identifying and managing trade-off i.e. favoring one property above the other. For instance, adding an undo function to an interactive system will improve usability by make it more efficient for users to recover from errors. However, adding undo functionality to a system increases significantly the number of lines of code and thus the likelihood of bugs. This paper focuses on dependability related issues and how dealing with them might bring additional concerns for the design of user interfaces and their associated interaction techniques. However, despite this specific focus on one property, similar constraints would apply to other conflicting properties.

Assuring that operators will be able to perform their activities even though the interactive system exhibits failures is one of the main issues to address when designing and implementing interactive systems in safety critical contexts. Exploiting methods, techniques and tools from the dependable computing field [START_REF] Laprie | Basic Concepts and Taxonomy of Dependable and Secure Computing[END_REF] can ensure this even though they have not been designed and developed to meet the challenges of interactive systems [START_REF] Fayollas | A Software-Implemented Fault-Tolerance Approach for Control and Display Systems in Avionics[END_REF]. Such approaches can be dived into two main categories:

-The zero-defect approaches (usually based on formal methods [START_REF] Weyers | The Handbook of Formal Methods in Human-Computer Interaction[END_REF]) that aim at guaranteeing that the interactive system will be defect free. While this has been proven a good mean for detecting and removing faults and bugs at development time, natural faults (such as bit-flips due to radiations) are beyond their reach.

-The fault-tolerant approaches that promote the use of redundancy (multiple versions of the system), diversity (the various versions are developed using different means, technologies and providers) and segregation (the various versions are integrated in the operational environment by independent means e.g. executed on different computers, using different communication means, …). Segregation ensures that a fault in one of the versions will not induce a fault in another versionusually called common point of failure.

One of the way to apply dependability principles to the user interface of the interactive system is to propose redundant user interfaces offering multiple ways for the user to perform operations. This can be displaying the same information on different screens or offering multiple input devices for triggering the same action. This can also be performed at the interaction technique level as presented in [START_REF] Navarre | A Formal Approach for User Interaction Reconfiguration of Safety Critical Interactive Systems[END_REF] where mouse failures were mitigated by the use of "similar" configurations based on use of multiples keys on the keyboard. However, to avoid user errors (such as capture errors [START_REF] Reason | Human Error[END_REF]) and increase learnability, it is important to ensure that the various user interfaces are "similar" at presentation and interaction levels. This concept of similarity has already been used in the field of web engineering [START_REF] Heil | Measuring and Ensuring Similarity of User Interfaces: The Impact of Web Layout[END_REF] but only with a focus of designing new web systems being consistent with legacy non-web systems. This paper refines the concept of similarity and shows how this concept is relevant at different levels of the architecture of interactive systems. The paper then presents a set of examples from the avionics domain where dependability is a major concern and where development of fault-tolerant mechanisms is a requirement from standardization authorities such as DO 178C standard [START_REF]C / ED-12C, Software Considerations in Airborne Systems and Equipment Certification[END_REF]. These examples present how similarity has been driving the design of multiple user interfaces even though they are as different as hardware only (interaction taking place through knobs and dials) and software mainly using WIMP interaction techniques. Conclusions and discussions for the workshop are presented in the last section.

Conflicts and Congruence between Similarity, Diversity and Redundancy in the area of interactive systems

In order to increase resilience to failures, fault-tolerance (i.e. guaranteeing the continuity of service), requires duplicated user interfaces for the command and control of a single system. This ends up with redundant user interfaces serving the same purpose. If those interfaces are built using the same processes and offer the same interaction techniques, it is possible that a single fault could trigger failures in both user interfaces. This could be the case for instance when using the idea of cloning the UI as proposed by [START_REF] Villanueva | Distributing web components in a display ecosystem using Proxywork[END_REF]. In order to avoid such common points of failure the redundant user interfaces must ensure diversity. Diversity can be guaranteed if the user interfaces have been developed using diverse means such as different programming languages, different notations for describing their specification, executed on top of different operating systems, exploiting different output and input devices, … Such diversity is only efficient if the command and control system offers confinement mechanisms avoiding cascading faults i.e. the failure of one user interface triggering a failure in the duplicated one. Such fault tolerant basic principles raise conflicting design issues when applied to user interfaces. Indeed, diversity requires the user interfaces to be very different in terms of structure, content and in terms of interaction techniques they offer, even though they must guarantee that they support the same tasks and the same goals of the operators [START_REF] Fayollas | Fault-Tolerant User Interfaces for Critical Systems: Duplication, Redundancy and Diversity as New Dimensions of Distributed User Interfaces[END_REF]. Another aspect is that they must be located in different places in the system i.e. distributed as this is one of the most efficient way of ensuring confinement of faults.

In that context, distribution of user interface does not concern the presentation of complementary information in different contexts (as presented in [START_REF] Martinie | A multi-formalism approach for model-based dynamic distribution of user interfaces of critical interactive systems[END_REF]) but the presentation of redundant information in those contexts.

In terms of design, it is important to be able to assess that the various user interfaces make it possible for the operators to reach their goals (this would be called similarity in terms of effectiveness). Beyond that, it is also important to be able to assess the relative complexity and diversity of these interfaces in order to be sure that operations will not be drastically degraded when a redundant user interface has to be used after a failure has occutred on another one. Studying the effective similarity (in terms of efficiency) at the level of input and output is thus required even though different type of displays and different types of input devices have to be used. This goes beyond the study of similarity at effectiveness level, but both contribute to the usability of the systems. It is important to note that all the other properties mentioned previously are intrinsic or extrinsic propeties of a given interactive system. Similarity is very special as it only has a meaning when two interactive systems are considered (or two different versions of a same interactive system). Such relative properties are usually less studied than absolute properties as the focus of interest is usually to favour a given property of a given system.

3

Examples from the Avionics Domain

The case study presents (in the area of aircraft cockpits) examples of redundant user interfaces. More precisely, we present in the context of the cockpit of the A380 (see Figure 1) aircraft. In this new generation of large civil aircrafts, the cockpit presents display units (that can be considered here as computers screens) of which some of them are offering interaction via a mouse and a keyboard by means of an integrated input device called KCCU (Keyboard Cursor Control Unit). Applications are allocated to the various display unit (DU).

In the A380, two redundant ways of using the autopilot are offered to the pilot in order to change the heading of the aircraft. One is performed using the electronic user interface of the Flight Control Unit (FCU on top of Figure 1)while the other one exploits the graphical user interface of the Flight Control Unit Backup interface and the KCCU (bottom of Figure 1).

KCCU FCU FCU Backup

Figure 1. Two possible means to control flight heading within the A380 interactive cockpit, one using the FCU and the other using the FCU Software application and the KCCU

Example one: entering a new value for heading

Figure 2 presents a zoomed view on the two means for entering a new heading of the aircraft. On the left-hand side of the figure, the editing of the heading is performed using a physical knob, which may be turned to set a heading value (this value ranges from 0 to 360). The selected value can be sent to the autopilot (called "engaged") by pressing the physical LOC push button below the knob. On the right-hand side, the heading is set using the keyboard of the KCCU and engaged by using the KCCU and its manipulator to click on the dedicated software LOC push button.

At a high level of abstraction (i.e. not taking into account the input and output devices), the task of setting a new value for the heading is the same on both user interfaces (they are similar at the effectiveness level). If described at a lower level, the description of these two tasks would be different, as they would require different physical movements from the pilots (they are thus not similar at the effectiveness level as for instance, the pilot would have to execute the FCUS application while the hardware FCU is directly reachable). It is important to note that there are other additional means to perform the same task (for instance controlling directly the aircraft using the sidestick) that are not presented here.

3.2

Example two: entering a set of parameters for the Navigation Display Figure 3 presents two different means to handle both barometer settings and parameters of the navigation display (NDpilot ND is the second screen on the left in Figure 1 while first officer ND is the second screen on the right). It illustrates how physical input devices (on the left-hand side of Figure 3) have been transposed into software components (right-hand side of Figure 3) handled using the KCCU (as in the FCUS presented in Figure 2). The general layout of both interface is quite close to that one, but the translation into a software application leads to different design options:  On the physical interface, the two barometer settings options (highlighted in yellow and on the bottom left part of both physical and software interfaces) are handled using two physical labelled push buttons (LS and VV) that are lighted on with a single light when the option is selected. The transposition of these two buttons in the software user interface results is a set of two software buttons that may be highlighted by changing the color of three horizontal lines. In this case, the two design options are quite similar.

 The General ND parameter settings (highlighted in green and on the top right part of both physical and software interfaces) are physically handled using physical push buttons without labels associated to labels displayed on a dedicated screen. These buttons behave in the same way as the two previous buttons. The software transposition is similar to the previous one, using both software push buttons and labels, and following the same layout constraints (relative position and size) as the physical interface.

 The Pressure editing (highlighted in red and located on the left-hand side of both physical and software interfaces) consists in the editing of a numeric value. The physical and software representations of this function follow two distinct design option. With the physical interface, this value is modified using a physical knob and the edited value is displayed on a dedicated screen while on the software transposition, this editing is performed using a classical text field that embed both editing and display of the value. It is thus possible on the software UI to use the arrow keys to navigate into the text box and modify one specific digit of the pressure, which is not feasible on the hardware UI.

 The ND range setting (highlighted in blue and on the bottom right part of both physical and software interfaces) is performed by selecting a range amongst a finite set of predefined values. In this case, the two design options are quite different too. On the physical interface, the task is performed using a knob that can rotate between the set of values, these values being physically written around the knob (making it visible at any time). The software translation of this interface is made up using a drop down combo box that embed both the display and selection of the value. In this case, the selectable values are only displayed while using the software component.

Example one: visualization of aircraft pitch and roll

Figure 4 presents two different design of the gyroscope instrument that aims at providing the pilot with information about the position of the aircraft relatively to the horizon (both pitch and roll). At the bottom right-hand side of Figure 4 the cockpit presents the physical analog display of these values. This device is also called the artificial horizon as the information it displays is similar to the view the pilots have when they look outside through the windshield. The software transposition of this instrument (on the left-hand side of Figure 4 called Primary Flight Display) embeds several other functions such as an altimeter or a speed controller. The graphical layout of the software UI is clearly inspired by the physical one which was, in the early days of aviation only a physical ball emerged in a container filled with liquid.

Research Directions

While the examples above focus on the presentation and interaction aspects of interactive systems, we are investigating other means to support similarity analysis and assessment also at user level:

-Investigating means of describing past experiences and practice of users to understand the level of familiarity a user may have with a given interaction -Investigating means of describing tasks (including knowledge, information and objects used to perform the tasks) such as with the HAMSTERS tool [START_REF] Martinie | Extending procedural task models by systematic explicit integration of objects, knowledge and information[END_REF] to assess similarity of tasks and goals -Investigating means of describing users' errors (including causes called genotypes and manifestation called phenotypes) in order to identify potential unexpected types of errors that could occur see [START_REF] Fahssi | Enhanced Task Modelling for Systematic Identification and Explicit Representation of Human Errors[END_REF].

In the area of aviation, the design driver for cockpit has been on targeting at similarity of command and displays even though this is not clearly stated. Indeed, looking at training, most airlines propose Cross Crew Qualification programs for pilots [START_REF]Cross Crew Qualification courses[END_REF]. As training is mainly based on tasks execution [START_REF] Martinie | Model-based training an approach supporting operability of critical interactive systems[END_REF] such a goals and task-based approach is critical and is the only way of designing and evaluating training programs evolutions.

Discussions and Conclusion

This paper has presented the similarity property for interactive systems offering redundant ways for the users to enter and perceive information. In order to ensure diversity and segregation (that are required for building dependable interactive systems) the similarity property may be violated. We have shown on the first example that the hardware and the software user interface are similar at the effectiveness level but distinct at interaction level. The following examples have shown bigger gaps in terms of similarity as the use of computing systems and graphical interfaces provides designers and developers with more advanced communication and interaction means. Digital devices are thus more informative and more efficient than the hardware ones. However, they are also less reliable than hardware systems and must not be used if failures are detected [START_REF] Fayollas | An Approach for Assessing the Impact of Dependability on Usability: Application to Interactive Cockpits[END_REF]. This means that the design and the evaluation of the training program is a complex and expensive activity requiring tools and technique to assess (and explain to trainees) gaps in similarity.

Figure 2 .

 2 Figure 2. Heading selection.

Figure 3 .

 3 Figure 3. Baro settings and Navigation Display configuration

Figure 4 .

 4 Figure 4. Physical and software representation of the aircraft gyroscope.