{'fraction_non_alphanumeric': 0.039774574617969, 'fraction_numerical': 0.014775477764531628, 'mean_word_length': 4.467410626111001, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Joint Motion Estimation and Source Identification using Convective Regularisation with an Application to the Analysis of Laser Nanoablations

We propose a variational method for joint motion estimation and source identification in one-dimensional image sequences. The problem is motivated by fluorescence microscopy data of laser nanoablations of cell membranes in live Drosophila embryos, which can be conveniently—and without loss of significant information—represented in space-time plots, so called kymographs. Based on mechanical models of tissue formation, we propose a variational formulation that is based on the nonhomogenous continuity equation and investigate the solution of this ill-posed inverse problem using convective regularisation. We show existence of a minimiser of the minimisation problem, derive the associated Euler–Lagrange equations, and numerically solve them using a finite element discretisation together with Newton’s method. Based on synthetic data, we demonstrate that source estimation can be crucial whenever signal variations can not be explained by advection alone. Furthermore, we perform an extensive evaluation and comparison of various models, including standard optical flow, based on manually annotated kymographs that measure velocities of visible features. Finally, we present results for data generated by a mechanical model of tissue formation and demonstrate that our approach reliably estimates both a velocity and a source.


Motivation
Motion estimation is a ubiquitous and fundamental problem in image analysis, see e.g. [5]. It is concerned with the efficient and accurate estimation of displacement fields in spatio-temporal data and has a wide range of applications, not necessarily limited to natural images. Optical flow [30] is one popular example of motion estimation, which designates the apparent motion of brightness patterns in a sequence of images and is based on the assumption of constant brightness. Recently, optical flow methods have been used for the quantitative analysis of biological image sequences on cellular and subcellular level. See, for instance, [2,9,10,21,31,34,39,44,50]. While the concept is in practice well-suited for natural scenes, the use of the less restrictive continuity equation, which arises from mass conservation, can be more favourable in certain scenarios. For instance, in [15,16] it is used for fluid flow estimation.
In developmental biology, the study of the morphogenesis of model organisms is specifically calling for image analysis methods that are able to extract time-dependent deformations and flow velocities from microscopy image sequences. Morphogenesis is the process that leads to an organism developing its shape as a result of the implementation of a genetic programme [28] and includes, among other mechanisms, tissue deformations. These deformations can be observed through video microscopy by fluorescently labelling molecules that are associated with compounds of mechanical relevance within an embryo [32].
In many cases, these molecules are organised spatially in discrete structures, such as cell membranes. To compute deformations on the level of these molecular structures, segmentation or detection, and subsequent tracking are the methods of choice [22]. As a result, detailed knowledge of the mechanics of morphogenetic processes can be gained [7]. In case the recorded tissue lacks structure, particle image velocimetry (PIV) is generally used to compute (sparse) displacement fields in image sequences. See, for example, [37,47].
One difficulty in abovementioned approaches is that the observed structures often have a short life time and are being degraded during the observation, while new structures of the same type are being created [47]. Indeed, these molecular structures, and thus their fluorescent signal response, can be described by an advection-reaction equation rather than pure advection [42]. The signal variations due to reaction are a source of error in the estimation of motion that we propose to address in this paper.
Beyond the need of measurements of dense velocities of moving fluorescentlylabelled molecular structures that prove robust with respect to the reaction term, it is of general interest to quantify the reaction term itself [42], and of general interest to extract information about all the physical quantities and processes, such as diffusion, that govern the observed tissue flow. For this, an accurate estimation of the velocity field corresponding to the time evolution of the distribution of fluorescent molecules is crucial [42,47]. The laser ablation is applied at frame number five, i.e. between the first and the second image shown. Its location is indicated with a magenta arrow. Observe the instantaneous recoiling and tissue loss of the cut membrane, and the subsequent growth of tissue in the cut region. Moreover, note the changes in contrast over time and the line artefacts.
In this article, we argue that utilising variational motion estimation can help to identify physical quantities by estimating velocities in real data. For simplicity, we choose to demonstrate this method in a case where the biophysical data can be reduced to one space dimension. This allows to create convenient space-time representations, so-called kymographs.
One-dimensional data are indeed relevant in tissue dynamics when the cellcell junctions are found to be aligned along a straight line called a supracellular actomyosin cable [8,40]. A common experiment to investigate the function of these cables is to cut them locally using intense laser illumination [25] and to observe the dynamics that follow. Figure 1 illustrates a prototypical two-dimensional (timelapse) fluorescence microscopy image sequence where a cable is being severed by such a laser ablation.
Most of the relevant dynamics occur along the cable itself, thus projecting the recorded signal in a narrow stripe of less than two micrometers along its average direction preserves most of the information. See Fig. 2 for an example of a kymograph obtained from the image sequence shown in Fig. 1. Variations in fluorescence intensity are clearly visible and displacements of features can easily be measured, e.g. with existing (tracking) tools [12,13,38,41]. See also Fig. 8 for an example of manually created tracks.  Fig. 1. A supracellular cable is indicated with a magenta rectangle. To reduce the problem to one dimension, the data is summed along the transverse direction within the rectangular zone. The right image depicts the kymograph obtained from this dimension reduction. Time runs from top to bottom and the horizontal black line at frame five indicates the time of the laser ablation, during which the signal acquisition was paused.
Analysing such data is challenging for many reasons. First, simultaneous estimation of velocity and decay or increase of the signal renders the problem ill-posed, as we will illustrate below, and suitable (qualitative) assumptions on favoured solutions are required. Second, the obtained kymographs are very noisy, contain artefacts due to the acquisition technique, and sometimes suffer from off-plane motion of the cables. Third, the velocity field potentially contains discontinuities at the time of the laser ablation. Fourth, data is missing during the application of the laser cut and only a limited field of view is available due to the nature of the kymograph. Moreover, bleaching of tissue leads to a decrease in contrast when being exposed over long periods of time. See Figs. 1 and 2 for illustration of these issues. In this work we address mainly issues one and two.
Motivated by the laser nanoablation problem we restrict ourselves to the onedimensional case and denote by Ω ⊂ R the spatial domain. For T > 0, we model the actomyosin concentration as a function f : (0,T) × Ω → R that is proportional to the observed fluorescence response. In the following, we assume that it solves the Cauchy problem for the nonhomogeneous continuity equation in one dimension: Here, v : (0,T) × Ω → R is a given velocity field, k : (0,T) × Ω → R a given source (or sink), and f 0 : Ω → R a given initial condition. Solution theory for problem (1) is closely related to solutions of the initial value problem via the method of characteristics, see e.g. [24,Chap. 3.2]. Here, the map φ : [0,T) × Ω → Ω denotes a so-called flow of v. Existence and uniqueness of problem (1) can be established by application of the Picard-Lindelöf theorem to (2), provided that both v and k are sufficiently regular. In particular, it can be shown that (2) has a unique (global) solution and, moreover, that this solution is a diffeomorphism. For an introduction and for details we refer, for instance, to [18,60].
In this work we are concerned with the solution of the inverse problem associated with (1), which is to estimate a pair (v, k) from (potentially noisy) observations f . Its ill-posedness is immediate because the problem is underdetermined and, given a for v 2 differentiable denotes a solution as well. In addition, it is easy to see that-without further assumptions on v and on k-the pair (0, ∂ t f ) is always a solution, albeit not a desired one.
In view of this ill-posedness and this ambiguity we consider the variational form where the data term D(v, k, f ) is the squared L 2 norm of the first equation in (1), R(v, k) a suitable regularisation functional, and α > 0 is a regularisation parameter. In this article we consider different choices of R. While source identification has been treated before in the literature, e.g. in [4], the main goal of this article is to recover a source (or sink) k which is constant along characteristics of the flow (2). In other words, we are interested in utilising the convective derivative for regularisation. Its use is inspired by the work in [33], where the convective derivative of the velocity along itself was used. From a physical perspective this choice seems natural as, in comparison to using, for instance, the L 2 norm or the H 1 seminorm for regularisation of k, it is consistent with the movement of the tracked cell tissue, which can be assumed to be the main origin of changes in the observed fluorescence intensity. However, from a numerical point of view this choice comes at the expense of having to solve nonlinear optimality conditions.

Contributions
The main contributions of this article are as follows. First, we propose a variational model based on the non-homogeneous continuity equation for joint motion estimation and source identification in kymographs. Second, we study the variational properties of utilising the convective derivative (4) for regularisation. Following [33], we establish existence of minimisers of the nonconvex functional. Third, for the numerical solution of the corresponding nonlinear Euler-Lagrange equations we propose to use Newton's method and a finite element discretisation. Fourth, we present numerical results based on kymographs of laser nanoablation experiments conducted in live Drosophila embryos. Moreover, we provide an extensive experimental evaluation of different data fidelity and regularisation functionals based on manually created tracks, and evaluate our approach using synthetic data generated by solving a mechanical model of tissue formation.

Related Work
In [30], Horn and Schunck were the first to pursue a variational approach for dense motion estimation between a pair of images. They considered a quadratic Tikhonovtype functional that relies on conservation of brightness and used (squared) H 1 Sobolev seminorm regularisation. This isotropic regularisation incorporates a preference for spatially regular vector fields. Well-posedness of the functional was proved in [51] and the problem was solved numerically with a finite element method. For a general introduction to variational optical flow see, for instance, [5].
In [57], the problem was treated on the space-time domain and extended to incorporate both spatial as well as temporal isotropic regularisation. In [56], a unifying framework for a family of convex functionals was established, and both isotropic and anisotropic variants were considered. We refer to [54] for more details on nonlinear diffusion filtering, and to [55] for an overview of numerous optical flow models and a taxonomy of isotropic and anisotropic regularisation functionals.
The convective derivative has already been used in several works. For instance, in [14] for simultaneous image inpainting and motion estimation. In [43], an optical flow term was incorporated in a Mumford-Shah-type functional for joint image denoising and edge detection in image sequences. Moreover, in [11] it was used for joint motion estimation and image reconstruction in a more general inverse problems setting. In [33] the convective acceleration was used for regularisation together with a contrast invariant Horn-Schunck-type functional. The corresponding nonlinear Euler-Lagrange equations were solved using a finite element method and alternating minimisation.
According to [16], the article [52] is credited for being the first to propose the use of the less restrictive continuity equation for motion estimation. Later it was used, for instance, to find 3D deformations in medical images [53,20], to analyse meteorological satellite images [6,16,61], and to estimate fluid [15,59] and blood flow [3] in image sequences. For a general survey on variational methods for fluid flow estimation see [29].
Whenever mass conservation is not satisfied exactly, e.g. due to illumination changes, it can be beneficial to account for these violations. For instance, in [27] they incorporated physical models. In [4] they simultaneously estimated image intensity, flux, and a potential source. In contrast to our work, only L 2 integrability of the source was assumed and the constraint was enforced exactly, leading to an optimal control formulation, which was solved with a finite element method.

Preliminaries
Notation For T > 0 and Ω ⊂ R a bounded, connected, and open set we denote by E = (0,T)×Ω the spatio-temporal domain and by ∂E its boundary. For a smooth function f : E → R we denote by ∂ t f , respectively, ∂ x f the partial derivatives with respect to time and with respect to space, and by ∇ f = (∂ t f , ∂ x f ) its spatio-temporal gradient. The space-time Laplacian of f is denoted by ∆ f . Analogously, for a smooth vector field w : E → R 2 with w = (w 1 , w 2 ) , its gradient is denoted by ∇w = (∇w 1 , ∇w 2 ) and its spatio-temporal divergence is given by ∇ · w = ∂ t w 1 + ∂ x w 2 . Moreover, we will write A B whenever there exists a constant c > 0 such that A ≤ cB holds. Finally, by the Cauchy-Schwarz inequality and application of Young's inequality, we have Here, and in the following, we use · L 2 instead of · L 2 (E,R) and · L 2 (E,R 2 ) for simplicity.

Convective Derivative
Let φ : E → Ω be a flow through the domain Ω, i.e. for every t ∈ (0,T) the map φ(t, ·) : Ω → Ω is a diffeomorphism and, for a fixed starting point x 0 ∈ Ω, the trajectory φ(·, x 0 ) is smooth. With every trajectory φ(·, x 0 ) that originates at x 0 ∈ Ω we can associate a velocity at every time t ∈ (0,T) via (2). Thus, a flow φ gives rise to a scalar (velocity) field v : E → R by means of where φ −1 (t, x) denotes the inverse of φ(t, ·) at x ∈ Ω, which is the starting point of the curve that passes through x at time t. For a scalar quantity k : E → R, we define the convective derivative of k along a flow φ, denoted by D v k, as For convenience we adopt the notationv = (1, v) . Clearly, D v k = ∇k ·v vanishes for pairs (v, k) such thatv ⊥ ∇k in R 2 . Moreover, from a straightforward calculation we obtain wherevv is the matrixvv As noted in [33], the action of D v k becomes more clear by writing (7) as The first identity states that, for fixedv, minimisation of the functional promotes functions k that vary little in the direction ofv. In addition, the weighting factor gives importance to regions of largev. On the other hand, the second identity states that, for fixed k, minimisation of (8) favours functionsv that are tangent to the level lines of k. This time, importance is given to regions where ∇k is large. See Fig. 3 for illustration. The connection to anisotropic diffusion [54] is immediate when considering the Euler-Lagrange equations associated with (8). They read Note that the system (9) is nonlinear.

Variational Model and Existence of a Minimiser
In this section we formulate the joint motion estimation and source identification problem and study the use of (8) as regularisation functional. For simplicity, we will denote a velocity-source pair by w = (v, k) . In the following, we intend to minimise a variational formulation of the form where R is yet to be defined and α > 0 is a regularisation parameter.
In order to fix a concrete choice of R let us discuss two issues. First, it should be selected such that the functional is well-defined for an appropriate function space. In particular, we require the weak derivative of v in the data term in (10) to exist and to be bounded with respect to an appropriate norm. Second, ∂ x ( f v) and k need to differ qualitatively in order to obtain a meaningful decomposition of the signal ∂ t f .
Before we state the model, let us follow the ideas in [33, Chap. 5.2.2] and discuss some issues arising with the choice To ensure well-definedness of the functional (10) we derive, by application of (5) and HÃ ¶lder's inequality, the estimates As a consequence, minimisation over the space seems appropriate for data f ∈ W 1,∞ (E) in the space of functions with essentially bounded weak derivatives up to first order. However, in order to show existence of minimisers of (10) by application of the direct method [19], one requires a coercivity estimate of the form where b ≥ 0 is some constant and the norm of the space X is defined via Without further restriction of the data f , the above functional is not coercive with respect to this space as the following example shows.

Example
Let f = ax for some 0 < a < +∞ and, consequentially, we have that ∂ x f = a. Let {w n } be the sequence with w n = (n, an) , for n ∈ N, and observe that w n X → ∞ as n → +∞, while the value of the functional E stays bounded, in fact, zero.
Moreover, it is clear that, for pairs w = (0, k) , the inequality (12) cannot hold in general. As a remedy, we consider minimising over all pairs w = (v, k) arising from the Sobolev space H 1 (E, R 2 ). Its norm is defined via In further consequence, our goal is to find a minimiser w ∈ H 1 (E, R 2 ) of the functional E in (10) with and regularisation parameters α, β > 0. Our final model thus reads Let us establish the existence of a minimiser to the problem min w ∈H 1 (E,R 2 ) E(w). The following lemma will be used to show coercivity of E and is along the lines of [51, p. 29]. holds.
Proof First, observe that, by the Cauchy-Schwarz inequality, the assumption ∂ x f const. is equivalent to where |E | denotes the measure of E.
Next, suppose to the contrary that there is no C > 0 such that (14) holds. Then, for all n ∈ N there exists w n ∈ L 2 (E, R 2 ) such that Letw n w n / w n L 2 withw n = (ṽ n ,k n ) and obtain ∇ −1 f ·w n together with (15) implies thatṽ n ,k n → 0 as n → +∞, contradicting the assumption. Therefore, (14) holds.
In particular, Lemma 1 holds for the (componentwise) average w E of w, defined as The following proposition is a straightforward adaptation of [33, Prop. 1] and utilises the direct method in the calculus of variations [19].
Proof The functional is proper and bounded from below since all terms are nonnegative and, for w identically zero, Next, we show coercivity of (13) with respect to H 1 (E, R 2 ). Observe that The chain of inequalities follows from (5), (14), the Poincaré-Wirtinger inequality [24,Chap. 5.8], and the assumption f ∈ W 1,∞ (E). Coercivity of E with respect to H 1 (E, R 2 ) then follows since ∇w 2 [24,Chap. 5.7]. As a consequence, there exists a subsequence, also denoted by w n , such that w n →ŵ in L 2 (E, R 2 ). Then, weak lowersemicontinuity of E follows by application of [19,Thm. 3.23] since, for fixed w, we have that all terms are convex in ∇w. In particular, |D v k | 2 is a quadratic form and therefore convex in ∇w, sincevv in (7) is symmetric positive semidefinite.
Let us add, however, that the convective regularisation functional F is nonconvex, as the following example shows. We pick w 1 = (0, x) and w 2 = (1, 0) , and obtain Therefore, E is nonconvex in general and several minima might exist. However, for β = 0 a unique minimiser exists.

Numerical Solution
In this section we derive necessary conditions for minimisers of (13) and discuss the numerical solution of a weak formulation by means of Newton's method.

Euler-Lagrange Equations
For convenience let us abbreviate where Id denotes the identity matrix of size two. Recall from Sec. 2.1 that ∆ and ∇ are spatio-temporal operators. Moreover, the natural boundary conditions at ∂E are given by where n ∈ R 2 is the outward unit normal to the space-time domain E.
Let us highlight two aspects of (16). First, note that the system is nonlinear in the unknown w = (v, k) due to the convective regularisation. Second, as already mentioned in Sec. 2.1, there is a connection of the second set of equations in (16) with anisotropic diffusion with the diffusion tensor given by αId+ βvv . The investigation of existence and regularity of solutions of (16) is left for future research.

Weak Formulation and Newton's Method
We minimise E by applying Newton's method to the weak formulation associated with (16) together with boundary conditions (17). It is derived as follows.
Multiplying with a test function ϕ = (ϕ 1 , ϕ 2 ) ∈ H 1 (E, R 2 ), integrating by parts under consideration of (17), and adding both equations leads to the following variational problem: with Here, V is the Gâteaux derivative ofvv at v in the direction δv and reads We solve the nonlinear problem (18) with Newton's method, which proceeds as follows. Starting from an initial solution w (0) (0, 0) we update the solution according to the rule where δw is the update and n ∈ N 0 . Computing δw in each step requires to solve a linear variational problem: with a(δw, ϕ) = DG(w (n) ; δw, ϕ), Note that the dependence on the current iterate w (n) is through G and DG, defined in (18) and (19), respectively. In our experiments we observed that (20) typically converges within only a few iterations.

Discretisation
We have implemented the weak formulation in FEniCS [1], which can also handle Newton's method automatically. The formulation (18) was discretised using multilinear finite elements. Since kymographs serve as input data we have discretised the rectangular space-time domain E with a triangular mesh based on the regular grid so that every vertex of the mesh corresponds to one pixel of the image f and to one pair of values of the unknown w. Moreover, in the implementation we penalise weak derivatives differently in space and time. This results in four regularisation parameters α i j , with i = {v, k} and j = {t, x}, for the H 1 seminorm in (13) and one additional parameter β for the convective regularisation. Due to the equivalence of norms the existence result in Prop. 1 still holds true.
Integrals are computed exactly with an appropriate Gauss quadrature, which is automatically selected by FEniCS. Similarly, since the image f is represented by a piecewise multilinear function, products of partial derivatives of f that appear on the right-hand side of (21) are automatically projected onto the correct space.
As termination criterion for Newton's method we used the default criteria of FEniCS with both the absolute and the relative residual of (18) set to 10 −10 . The maximum number of iterations was set to 15. Convergence was typically achieved within only a few iterations, which usually amounted to just a few seconds of computing time on a standard consumer laptop. It needs to be mentioned that, in our experiments we found that, the method fails to converge when the parameters α k t and α k x are chosen too small in comparison to β. This is, however, in line with the theoretical results in Sec. 2.2, which require H 1 regularity of k.
Both the source code1 of our Python implementation and the microscopy data2 used in the experiments are available online.

Experimental Results
In this section we present numerical results. In the first part, we demonstrate the importance of estimating a source based on synthetic data. Then, in the second part we show results for nonsynthetic microscopy data. After briefly discussing the data, we investigate the effects of varying the regularisation parameter β, which controls the amount of convective regularisation. We then qualitatively compare the (standard) 1D variational optical flow model with our continuity equation-based formulation for several choices of the regularisation functional αR in (10). In addition, we present a quantitative evaluation of the considered models based on recoil velocities obtained from manual tracking of features in kymographs. Finally, in the last part, we evaluate our approach based on data coming from the solution of a mechanical model of tissue formation.
In all results, the computed velocity fields are presented visually with the help of streamlines, see e.g. [58]. These are integral curves computed by numerically solving the ordinary differential equation (2) for the estimated velocity v and a selected number of initial points. The resulting curves are then colour coded according to their velocities and shown superimposed with the kymograph data. This representation is more comprehensible and allows to visually check whether the estimated velocities are approximately correct.

Analytical example
In order to demonstrate the necessity of estimating a source when the changes in the signal cannot be explained using mass conservation we conducted experiments for synthetic data. To this end, we generated a signal f , given as on a periodic domain Ω = (0, 1) and for the time interval [0, 1]. The parameters were set to τ = 1 and to λ = 1/(4π). This signal shifts to the right with constant velocity v 0 = 0.1 and decays exponentially over time in its magnitude. It can easily be verified that the source is given by k(t, x) = − f (t, x)/τ. We then solved the variational problem using in one case the homogenous and in one case the non-homogenous continuity equation. Periodic boundary conditions in space were enforced, and regularisation parameters were set to α v j = 10 −3 , α k j = 10 −4 , for j ∈ {t, x}, and to β = 10 −3 . Since k is not constant along characteristics in this example we haven chosen β quite small. Figure 4 illustrates the results of this experiment. It can clearly be seen that not accounting for the decay of the signal leads to a velocity that is significantly different to v 0 in this example, whereas using the non-homogenous continuity equation estimates both the velocity and the source (not shown) very well. Fig. 4: In this example, we demonstrate the necessity of source estimation for decaying data. Here, we have generated a signal that shifts to the right with constant velocity v 0 = 0.1 and simultaneously decays exponentially (in magnitude). Shown are the signal f and streamlines computed from the estimated velocities using the homogenous (left) and the non-homogenous (right) continuity equation.

Microscopy Data and Acquisition of Kymographs
The data at hand are 2D image sequences of living Drosophila embryos recorded with confocal laser-scanning microscopy. We refer to [48] for the used microscopy technique and for the preparation of flies, as well as for the details of the laser ablation method. For this study we recorded 15 image sequences, all of which feature cell membranes that have been fluorescently labelled with E-cadherin:GFP, see [48]. The image sequences feature a square region of approximately 42.2 × 42.2 µm 2 at a spatial resolution of 250 × 250 pixels. A typical sequence contains between 60 and 100 frames that were recorded at a temporal interval of 727.67 ms, and the recorded image intensities f 2D+T are in the range {0, . . . , 255}.
Each of the sequence shows a single plasma-induced laser nanoablation, which led to the controlled destruction of tissue in a linear region of 2 µm length that is approximately orthogonal to the cable. This ablation is expected to have a width of the order of the size of one pixel. Recall that in Fig. 1 we show such a typical dataset.
In order to obtain a kymograph from each microscopy sequence, we first labelled the location of the intersection between the ablation line and the actomyosin cable with a point c ∈ R 2 . Then, we visually determined an approximate orientation, given by a unit vector e ∈ R 2 , of the selected cable by defining a straight line of length 2L + 1 pixels which passes through c. Typically, L = 100 pixels is sufficient for the considered datasets.
To create a one-dimensional image sequence we used standard nearest neighbour sampling along the abovementioned line. A nearest neighbour in terms of the pixel locations P of f 2D+T is given by Then, the sampling points that are separated by distance one and lie along the abovementioned straight line are given by p i = c + ie, for i = −L, . . . , L. A nearest neighbour interpolation of f 2D+T along this line is given by f 2D+T (t, N(p i )).
However, due to noise, the spatial extent, and minor displacements of the cable in orthogonal direction, we also considered sampling points that lie on the parallel line of distance j = −h, . . . , h. These points are given by Then, at i = 1, . . . , 2L + 1, we define the intensity of a kymograph f δ (t, i) as In other words, the intensity at i is given as the sum of nearest neighbour interpolations at points that lie on an orthogonal straight line. The superscript δ indicates noise in the created kymograph. We found that h = 5 pixels generates satisfactory kymographs.
For this process, we used the reslice tool in Fiji [49] with a slice count of 2h and no interpolation selected, after manually placing a straight segment along a selected actomyosin cable, see Fig. 2 for illustration. Subsequently, a projection with the SUM option selected creates the final kymograph.
Since the image acquisition is paused during the laser ablation, see Fig. 2, we simply replaced the missing frame with the previous one. Moreover, we applied a Gaussian filter to the kymograph f δ to guarantee the requirements specified in Sec. 2.2 and scaled the image intensities to the interval [0, 1]. The kernel size of the Gaussian filter was chosen as 10 × 10 pixels and the standard deviation set to σ = 1. The filtered and normalised kymograph is denoted in the following by f .

Qualitative Comparison
In the first experiment, we investigated the effect of the convective regularisation for one chosen kymograph. To this end, we solved the necessary conditions (16) as outlined in Sec. 3 for varying regularisation parameter β. Figure 5 shows streamlines for the estimated velocities together with the computed sources. Since we did not find any significant difference between the results for β = 10 −4 and for β = 0, we simply omit the latter.
The two main findings of this experiment are as follows. First, for the chosen dataset, the convective regularisation can help to estimate more accurate velocities shortly after the ablation was applied. As can be seen in Fig. 5c (left) and Fig. 5d (right), this leads to a more accurate estimation of the recoil velocities at the cut ends. In Fig. 5d we display a magnified view of the results in the cut region.
Second, as expected, with increasing β the estimated source gets more regular in the direction of the flow, see Figs. 5a-5c (right). In particular, the oscillations in space and time in the estimated source, which can be most likely attributed to noise and to artefacts created during the acquisition, decrease significantly. This certainly can lead to better interpretability and help to get a better understanding of the estimated reaction term k. Observe also the decrease in the magnitude of k as β increases. This is made apparent by the colour coding.  In the next experiment, we qualitatively compare the standard variational optical flow model with H 1 seminorm regularisation and the continuity equation-based model in (10) paired with different regularisation functionals. The first model we consider is based on the optical flow equation [30] in one space dimension, which reads and assumes that f is constant along characteristics of the flow (2). Even though no source k is estimated in this model, the main motivation for its inclusion in the evaluation is that it can serve as a baseline method. The quantitative evaluation in Sec. 4.4 is based on manually created tracks that follow highly visible features in (a) Streamline representation of the velocity obtained by solving the variational optical flow problem (23) with spatio-temporal H 1 seminorm regularisation. Parameters were set to α v j = 5 · 10 −3 , for j ∈ {t, x }.
(d) Result obtained using αR CR with parameters set to α v j = 5 ·10 −3 , α k j = 10 −4 , and β = 2.5 ·10 −3 . Fig. 6: Qualitative comparison of different models based on one chosen dataset. Regularisation parameters were chosen manually so that the recovered streamlines best matched the cut ends.
the kymographs and, in many cases, roughly preserve their intensity as the tissue deforms. See Fig. 8 for illustration. We highlight that, given ∂ x f 0, equation (22) admits a unique solution, namely v = −∂ t f /∂ x f . However, due to noise degradation and aforementioned artefacts it is beneficial to solve (22) in a variational framework. Therefore, we consider minimising the functional analogously to the solution method outlined in Sec. 3. In contrast to the functional (13), the corresponding Euler-Lagrange equations are linear and the weak formulation can be solved directly. All other functionals we investigate are based on (10) and read Here, αR i represent different choices in the regularisation and, consequentially, also in the function space we minimise over. Since the above functional doesn't require any Sobolev regularity of the source k, we investigate the setting where v ∈ H 1 (E) and k ∈ L 2 (E), that is with γ > 0. Moreover, we also consider the choice v ∈ H 1 (E) and k ∈ H 1 (E). The regularisation functional then reads Finally, we will present results for the setting v ∈ H 1 (E) and k ∈ H 1 (E) together with the convective regularisation. This is the main model investigated in Sec. 2.2 and is given by In Fig. 6 we show minimising functions for models (23)- (26) and in Fig. 7 we display a magnified view of the cut region. For all models, the regularisation parameters were chosen manually so that the streamlines obtained from the computed velocities best matched the cut ends which resulted from the laser ablation. For comparison, we also refer the reader to Fig. 8, which shows manually tracked cut ends. The four main observations from this experiment were as follows.
First, in Fig. 6a, which shows a minimising function of the optical flow model (23), it is clearly visible that many characteristics follow paths of constant fluorescence intensity. See also the two outermost markers in Fig. 7 (left). However, this model underestimate the recoil velocity shortly after the cut and may lead to wrong characteristics, see the middle marker in Fig. 7 (left).
Second, in Fig. 6b we illustrate a minimising pair (v, k) of the continuity equation-based model using αR H 1 −L 2 . It is apparent that k captures both noise and artefacts, and leads to an undesirable underestimation of the velocities in the cut region. See also Fig. 7 (second from left). While this model unsurprisingly results in the smallest residual error (in norm) of the non-homogenous continuity equation, cf. also Table 1, it is insufficient for a meaningful quantification of tissue loss and growth.
Third, in Fig. 6c we depict a minimising pair for the model using αR H 1 as regularisation functional. Less noise is picked up by the source k and the recovered velocities v are closer to what one would expect in the cut region, see also Fig. 7 (third image from left). However, undesired oscillatory patterns are present in the source k.
Finally, in Fig. 6d we show a minimiser of the model using convective regularisation, that is αR CR . This particular choice leads to significant visual improvement both in the recovered velocity and the estimated source. Specifically, as the streamlines in the cut region right after the laser ablation in Fig. 7 show (see top markers), the movement of the tissue is captured very well. In addition, bright features such as the left cut end, cf. the bottom marker in Fig. 7, are followed accurately. Moreover, the changes in the fluorescence intensities are indicated nicely in the visualisation of the source k, see Fig. 6d (right). In particular, the significant increase between the cut ends towards the end of the sequence, possibly due to wound healing, is indicated adequately.
In summary, it can be said that our variational model based on the nonhomogenous continuity equation together with convective regularisation can lead to improved results compared to existing models when parameters are selected by hand and comparison is performed visually.

Quantitative Comparison based on Measured Recoil Velocities
In this section, we compare the models (23)- (26) presented in the previous section from a quantitative point of view. Our evaluation is based on manually measured recoil velocities of clearly visible features in the kymographs, and include the cut ends resulting from the laser nanoablation.
For this comparison, we have annotated 15 kymographs in Fiji [49] to obtain discrete trajectories to compare to. See Fig. 8 (left) for an example. These tracks can then be used to compare to either the estimated velocities or to the computed characteristics that solve (2) numerically. Observe in Fig. 8 that all trajectories start only after the ablation, which is the main region of interest from a tissue mechanics point of view.
Before presenting the comparison, let us briefly discuss the methodology and the used evaluation criteria. Since some created tracks do not feature a coordinate for each time instant we interpolated each track i with a third-order spline φi. As a result, velocities ∂ t φ i can be computed conveniently and used for comparison to the velocities estimated by the variational approach. See Fig. 8 for an example of tracks (left) and their corresponding spline interpolations (right), which are colour coded according to their velocities.
In our experiments we found that each kymograph requires the regularisation parameters to be adjusted individually. Therefore, in the experimental comparison, we performed a search over all parameter combinations for whichever parameters are applicable to the respective model. In the case of the most complex model, i.e. the one that includes the convective regularisation (26), this amounted to probing 5 5 parameter combinations per dataset. For each of the criteria listed below, we recorded the best result that was obtained with each model and for each kymograph during the parameter search.

Error in Residual
In the first comparison, the goal was to see which model best fits the recorded data.
In Table 1 we report the L 1 norm of the smallest observed residual of the underlying model equation, that is As expected, the model αR H 1 -L 2 results on average and for every single dataset in the smallest residual. However, as already mentioned in Sec. 4.3, due to the high noise level it is only of limited use for quantifying the reaction term as k captures a significant amount of noise. The main finding of this experiment is that the optical flow model is by far not capturing the entire essence of the dataset, which is indicated by the high residual in comparison to the continuity equation-based models. Moreover, let us also highlight that the residual is on average not significantly increased in comparison to the αR H 1 model when the convective regularisation is used in addition.  Table 1: Residual error for every dataset and every model investigated. For each dataset we show the smallest residual in L 1 norm that is obtained during the parameter search. Fig. 9: Plotted is the mean error (28) (vertical axis) for the model αR H 1 for each of the 15 datasets (in different colours) and all probed parameter combinations (horizontal axis). It can be seen that no single tested parameter settings led to a small error for all datasets.

Error in Velocity
In the second comparison, we looked at the absolute error between the velocity of each manually created track and the velocities estimated with our models. For a particular track φ i of a dataset, we define this error at time t ∈ [0,T i ] as Here, T i > 0 is the length of the track and we have assumed for simplicity that all trajectories start at t = 0. The velocity v needs to be interpolated, since φ i is a spline representation.
In further consequence we computed, for each dataset and for each parameter configuration, the mean squared L 2 norm of the error in velocity along all its N tracks. It is given by In our experiments, in none of the tested models we could find a single parameter combination that worked well for most datasets in terms of the mean error E(v). This can be seen, for example, in Fig. 9, where we have plotted exemplary for each dataset and for each parameter setting in (27) the error E(v) for the model αR H 1 .
As a consequence, we report in Table 2 (left) the best mean error E(v) that we obtained for each dataset by the grid search. The main findings are as follows. First, and most importantly, the continuity equation-based model with convective regularisation performed on average as well as the optical flow-based model when using E(v) as evaluation criterion, with the advantage of simultaneously yielding an estimate of the source.
A possible explanation for the comparably good performance of the optical flowbased model is that the manually created tracks approximately constitute trajectories of constant intensities rather than the characteristics associated with (1). However, in combination with the findings presented in Table 1, which show that the average residual is much smaller when using a continuity equation-based model with H 1 seminorm or convective regularisation, we are confident to state that these models are capable of estimating a meaningful source that can explain significantly more details of the observed signal.
In addition, we also evaluated (28) with · 2 L 2 replaced by · L ∞ , see Table 2 (right). Qualitatively, this leads to slightly different results for some datasets but still supports our main findings.
Let us illustrate the advantage of the convective regularisation on the basis of one particular dataset. In Fig. 10 we show the best result obtained for dataset number eight for three models. For this particular dataset, the model αR CR outperforms all other models according to Table 2. Figure 10a (left) shows the best result for the optical flow model and Fig. 10a (28), and the maximum error (right) for every dataset and every model investigated. For each dataset we show the smallest error that is obtained during the parameter search.
velocities. Notice the inaccurate velocity between the cut ends shortly after the laser ablation. Moreover, towards the end of the sequence, where the both cut ends meet again, the characteristics seem inappropriate. In Fig. 10b we display the result for the model αR H 1 . As can be seen in Fig. 10b  (left), the estimated velocity is improved significantly in the cut region and in the problematic region towards the end of the sequence.
In Fig. 10c (left) we show the result for the model including the convective regularisation. i.e. for αR CR . The estimated velocity appears visually as good as in the previous model with the additional advantage that it allows a larger magnitude shortly after the laser ablation. Observe that in both cases the estimated source is both spatially and temporally very regular, and apart from beta the same parameters α i j were selected. However, in Fig. 10c (right) the effect of the anisotropic regularisation is clearly visible in the cut region.
Moreover, in Fig. 11 we show the best obtained results for two other datasets. For comparison, the top row shows the same dataset as in Figs. 5 and 6.
Finally, let us mention that, for the model αR CR , which we solve via Newton's method, in total only 53 computations did not converge out of the 3125 parameter combinations tested on 15 different datasets.

Comparison based on a Mechanical Model of Tissue Formation
In Sec. 1 we have motivated the use of the non-homogenous continuity equation (1) mainly through mechanical models that are known to describe tissue formation, (a) Streamline representation of the velocity obtained by solving the variational optical flow problem (23) with spatio-temporal H 1 seminorm regularisation with parameters α v t = 10 −1 and α v x = 5 · 10 −2 .
(b) Result obtained using αR H 1 with parameters α v t = 10 −1 , α v x = 5 · 10 −2 , α k t = 10 −1 , and Fig. 10: Best results obtained for dataset number eight in terms of the mean error (28) in the estimated velocity, see Table 2 (left). for example, in Drosophila. In order to see whether our variational formulation can reliably estimate velocity and source that both stem from such a process, we have implemented the partial differential equation-based model proposed in [26] to create synthetic data. In this model, the triplet (m, v, σ) solves the system for (t, x) ∈ (0,T) × (0, 1), subject to the initial and boundary conditions We refer to [45] for a derivation of this model. Briefly, (29) is an advection-reaction equation modelling mass conservation of myosin molecules m, subject to rates of adsorption k on > 0 and desorbtion k off > 0, Fig. 11: Best results obtained for dataset number three and five in terms of the mean error (28) in the estimated velocity, see Table 2 (left).
which we assume to be constants. The advection is determined by a mechanical problem, where mechanical balance (30) involves the stress σ(t, x) in the actin at the cell-cell junction and a friction force ξv exerted by the surrounding material, which we assume to be proportional to the velocity. Here, ξ > 0 is the coefficient of viscous drag. Finally, a constitutive model (31) of the junctional actin is proposed following e.g. [35,23,26] and involves viscous stresses and a pre-stress χm generated by myosin molecules. Again, η, χ > 0 are constants. In addition, (32) enforces zero flux at the spatial boundaries, and m 0 is an initial concentration. In our experiments we set it to While such a model captures the essential features of actomyosin behaviour [35,23,47], its mean field approach means that away from the perturbation caused by the laser cut, the concentration m will equilibrate to the trivial solution (m, v, σ) = (k on /k off , 0, χk on /k off ), which means that these parameters are uniform in space, and there will be no feature to track for an image analysis technique. This is also in contradiction with the experimental observations, where some material points along the cell-cell junctions exhibit accumulations of myosin that persist over time. One possible biophysical explanation for these accumulations is a locally larger density of actin binding sites.
This can be incorporated in the above model by introducing an additional variable ρ(t, x) and a constant k 0 off > 0 that modulates the off rate of myosin In addition, the density ρ obeys the conservation equation and satisfies an initial condition. In our experiments we set the initial ρ at t = 0 as ρ 0 (x) 1 + 1 + sin(40x + cos(40x)) 10 .
We solve the system (29)-(31) numerically, under additional consideration of (29), with a standard upwind finite volume discretisation paired with the forward Euler method. See, for instance, [46, Appx. B] for a brief description. For completeness, we briefly outline our implementation here.
We discretise the space-time domain [0,T]×[0, 1] using N t and N x equally spaced discretisation points in time and in space, respectively. For the discretisation of the unknowns we make use of a centred and a staggered grid, denoted by G c and G s , respectively. They are defined as where ∆t = T/N t and ∆x = 1/N x . The concentrations and the stress are then discretised on the centred grid and the velocity on the staggered grid, leading to In further consequence, the finite volume discretisation, see e.g. [36,Chap. 4], of (29) reads where F(m)| j± 1 2 is the flux at the cell boundaries, that is, at the nodes of the staggered grid. Assuming that m is constant on each cell and using the upstream value of m, the flux F(m)(t)| j± 1 2 F (i) j± 1 2 at time t = i∆t at the boundaries can be written as Similarly, we obtain F (i) j− 1 2 and, moreover, (17) results in zero flux at the boundaries.
Then, approximating d/dt with forward finite differences yields Analogously, an update equation for ρ is obtained.
The velocity at the current time step (i) is determined as follows. Using (30) to obtain v = ∂ x σ/ξ and substituting it into (31) yields the second-order elliptic partial differential equation for the stress σ. Here, the zero Neumann boundary conditions follow from (32).
In each time step we solve (36), given the concentration m (i) from the previous time instant, with a standard finite-difference scheme using centred differences on G c . Then, with the help of (30) the velocity at nodes j + 1/2 can be approximated with where we also use the fact that v is zero outside the spatial domain. Finally, the concentration m (t+1) i is updated according to (35) and ρ to its corresponding equation, and the procedure is repeated for the next time instant. The time interval is adjusted in each step so that the Courant-Friedrichs-Lewy condition is satisfied. This typically leads to intermediate results that are not recorded.
In our experiments we set N t = 300 and N x = 300, and the parameters controlling the time stepping were set to T = 0.1 and to ∆t = 2.5 · 10 −6 . The mechanical parameters in (30) and (31) were chosen as η = 1, ξ = 0.1, and as χ = 1. The parameters in (31) related to the source were set to k on = 200 and to k 0 off = 10. In order to evaluate our approach described in Sec. 2.2, we conducted two experiments based on this mechanical model. In the first experiment, we used the solution method outlined above to generate a concentration m, which was then used as input to our variational formulation defined in (13). In the second experiment, we generated a concentration by solving (29) for a set velocity v, effectively removing the mechanical part of the model.

Unknown Velocities
In this experiment, we solved the mechanical model (29)-(31) together with (33) numerically as outlined above. In order to simulate a laser ablation, the concentration m 0 is set to zero at nodes within the interval [0.495, 0.505]. In this way, a disruption (or loss) of concentration is simulated. Figure 12 (top) shows the solution (m, v) and the resulting source k = k on − k 0 off ρm. We then solved (13) with αR CR numerically based on the generated concentration. This is achieved by setting f m. However, in order to match the boundary conditions in (32) we also used zero Dirichlet boundary conditions for v in (18) at x ∈ {0, 1} and at t = 0. In Fig. 12 (bottom), we depict an approximate minimiser. The parameters were set to a v j = 10 −4 , a k j = 5 · 10 −5 , and to β = 10 −6 . Observe that both the velocity and the source are estimated approximately and are within the correct order of magnitude. However, let us add that the estimated source appears quite regular in comparison to the simulation, even though the regularisation parameters α k j were set comparably small.  (10) with αR CR . The left column shows the myosin concentration m together with streamlines obtained from the velocities, the middle column depicts the velocity v, and the right column illustrates the source k.

Predefined Velocities
In the next experiment, we removed the mechanical part and solved just (29) and (33) with a predefined velocity field v that could potentially resemble a laser nanoablation in cell membranes as pictured, for example, in Fig. 2.
We generated a velocity profile as follows. First, we define a characteristic φ : [0,T] → R that is supposed to follow a cut end via the ordinary differential equation where the constant v 0 > 0 is the initial velocity at time zero. Integrating with respect to time yields the integral curve where c 0 ≥ 0 is a constant defining the starting point of the curve. We then define a velocity fieldv so thatv(t, ·) is linear in the interval [0, φ(t)] and decays exponentially in (φ(t), +∞).
Here, τ, > 0 are constants and control the decay. Finally, we shift the origin to 1/2,  (29) and (33) given a set velocity v. The bottom images depict an approximate minimiser w = (v, k) obtained using (10) with αR CR . The left column shows the myosin concentration m together with streamlines obtained from the velocities, the middle column depicts the velocity v, and the right column illustrates the source k.

reflectv, and obtain
In Fig. 13 (top) we illustrate the solution (m, v, k) for this scenario, where v is set as in (37). The parameter c 0 was set to c 0 = 0.05 to match the width of the simulated laser ablation. The other parameters were set to v 0 = 1, τ = 0.075, and to = 0.05. All other settings and parameters were kept as in the previous experiment.
Then, we minimised (13) together with convective regularisation numerically based on the resulting concentration m. Figure 13 (bottom) shows the estimated velocity and source pair for the synthetic data. A similar behaviour as in the previous experiment can be observed.

Conclusions
In this article, we have investigated a variational model for joint velocity estimation and source identification in challenging fluorescence microscopy data of live Drosophila embryos that show the controlled destruction of tissue. We exploited the fact that a large proportion of tissue deformation occurs along one space dimension and allows to create kymographs. Our formulation is grounded on one-dimensional mechanical models of tissue formation and is based on the nonhomogenous continuity equation. We have discussed the ill-posedness of this problem and devised a well-posed variational formulation using convective regularisation of the source. Moreover, we have shown the connection of convective regularisation of the source to anisotropic diffusion. In a thorough experimental evaluation, we have demonstrated that motion estimation can benefit from simultaneously estimating a source and that convective regularisation may help to estimate velocities more accurately. Our numerical results show that this method could potentially help to quantify the reaction term in biological models of tissue formation. The extension of our models to more than one space dimension is left for future research.