Damien Besancenot
email: damien.besancenot@parisdescartes.fr

Radu Vranceanu
email: vranceanu@essec.fr.

AN EQUILIBRIUM SEARCH MODEL OF THE FRENCH DUAL MARKET FOR MEDICAL SERVICES RESEARCH CENTER

Keywords: JEL Classi…cation: D83, I11, I18 Equilibrium search, Medical fee dispersion, Dual market, Regulation

The French market for specialist physician care has a dual structure, including a sector 1 with regulated fees, and a sector 2 where physicians can freely choose fees. Patients who undergo a sequential search process for the best medical o¤er develop a reservation fee decision rule. We analyzed physicians'decisions to work in sector 1 or in sector 2, and their choice of fee in sector 2. The model features several pure strategy equilibria that can be ordered with respect to patient welfare. Policy implications follow.

Introduction

Out of the 223,571 physicians active in France in 2016 as recorded by the Statistical O¢ ce of the Health Ministry (DRESS), 47,4% were self-employed, and another 10,8% operated under a mix of self-employment and contract based services. 1 There is no ex-ante allocation of patients to physicians; patients are free to choose their physician. In particular, in the market for specialists, patients will search for the physician that best matches their expectations. A similar search process can be observed with respect to newcomers in a city who search for a general practitioner. This paper aims to analyze how trade frictions can alter the equilibrium fees posted by specialist physicians in a decentralized market for medical services, and to identify the implications for patient welfare, given the speci…c organization of the French health care sector.

Since 1980, self-employed specialist physicians in France have been allowed to choose between performing their medical practice in a heavily regulated "sector 1" or under the much looser regulation in "sector 2". 2 In 2014, 43% of the specialist physicians were working in sector 2, and this proportion is steadily increasing. 3 Fees charged by physicians in sector 1 are determined by the National Health Insurance in an agreement with physicians, and they are covered in full.

The …rst part of the fee is covered by the universal public insurance (better known as the "Social Security") up to a predetermined amount. The remaining amount is covered by complementary insurance as provided by several private insurance companies, to which more than 96% of the population had access in 2010. [START_REF] Besancenot | Quality and price dispersion in an equilibrium search model[END_REF] Physicians working in sector 2 can freely choose their medical fees. According to Article 53 of the Public Health Law, they should set their fees with "a sense of measure". Nonetheless, fees charged in sector 2 are increasing every year; in 2014, the premium charged in sector 2 compared to sector 1 represented on average of 23 euros/medical care, which 1 Data extracted on Mai 26, 2017 on http://drees.social-sante.gouv.fr.

2 Between 1980 and early 1990, general practitioners could also opt for working in sector 2. Carrère (1991) reveals that 80% of the GPs decided to remain in sector 1. As a consequence of these past choices, 6% of the total GPs were still working in the sector 2 in 2012.

3 See CISS Collectif Interassociatif sur la Santé, Observatoire citoyen des restes à charge en santé / 21 mai 2015, www.santeclair.fr/web/site/santeclair.fr /…les/ 150521_dp_ depassementshonorairesmedicaux -toujourshausse_1.pdf is 35,7% of the total fee. In general, patients will pay for the di¤erence between the fee in sector 2 and the fee in sector 1 from their own pockets, although private insurance can also cover part of it (Clerc et al., 2012).

Moving from sector 1 to sector 2 is an important decision for the physician. In addition to the possible negative impact on the number of forthcoming patients, moving to sector 2 comes with additional costs. Access to sector 2 is restricted to physicians with a proven high skill level, as attested by experience or prestigious degrees. Furthermore, physicians in sector 1 receive a subsidy for their contribution to the retirement system which is not granted to physicians in sector 2. Finally, sector 2 physicians will in general bear some higher "representation costs", such as a better location or a more expensive o¢ ce.

As more and more physicians chose to move to sector 2, patient associations are complaining that …nding a doctor in sector 1 has become extremely hard, and that the excessive fees charged in sector 2 make access to medical care increasingly di¢ cult. 5 These worrisome trends prompted in 2013 the Health Ministry and the National Health Insurance to encourage physicians in sector 2 to sign special agreements (referred to as Contrats d'Accès aux Soins), where they voluntarily agree to cap their fees (on average) and in turn receive additional bene…ts. These contracts further evolved in 2017, and are now referred to as OPTAM (Option Pratique Tarifaire Maitrisée) agreements.6 While this scheme is relatively complex, in brief, physicians can agree to charge fees that are at most twice the fees in sector 1 for a predetermined percent of their service. In turn, private complementary insurance is allowed to cover these "moderate" fees in full. This paper uses an equilibrium search model to analyze the choice of sector and price decisions of self-employed physicians, and the consequences of these choices for patients. The use of a search model appears as a natural choice, as the information acquisition costs are quite high in the market for medical services, in particular for those services that involve complex diagnostic and treatment procedures (Rochaix, 1989; Mehrez et al., 1995; Emons, 2001) 7 .

Patients are assumed to be heterogenous with respect to their willingness to pay for the medical service and search for the best medical o¤er in a sequential mode. Physicians decide on the sector of activity. Furthermore, should they chose sector 2, they must decide on the best fee level.

The key policy variable is the regulated fee in sector 1, which traditionally has been a subject of controversy between the patients'association, the medical unions, the National Health Insurance and the government. The model presents several pure strategy equilibria that can be ranked with respect to patient welfare. In many cases, physicians can extract a signi…cant part of the patients' surplus. Dunn and Shapiro (2014) revealed that in the US, medical fees are higher in areas where the spatial concentration of physicians'is higher. Our model shows that trade frictions can increase the market power of physicians even if the market concentration is low. This model will allow us to shed some light on the consequences of imposing caps on the sector 2 fees as recently implemented in France.

As a limitation of this short paper, the quality of the medical service is assumed to be invariant.

An extension which would allow for a heterogenous health care quality could be developed along the lines of the equilibrium search model with price and quality di¤erentiation by Besancenot and Vranceanu (2003). In addition, in building the model we assumed that the cost of moving from sector 1 to sector 2 is identical for all physicians. If physicians have di¤erent costs form moving from sector 1 to sector 2 the model could generate a mixed strategy equilibrium where some physicians choose the sector 1 and the others chose the sector 2. However, in this more realistic setting the essential decision on whether physicians in sector 2 post high or low fees would be guided by the same rationale as in the simpli…ed setting.

The model 2.1 Main assumptions

We consider the market for a given medical service. The problem features a unit mass of patients and a unit mass of physicians. In line with classical assumptions in the literature, patients engage in a sequential search process, meeting one physician per period (McCall 1971). We assume that search costs are merely the psychological cost of delaying the cure process, and involve a discount factor < 1. In the context of medical services, this discount factor incorporates the probability that the disease aggravates quickly and can thus be very low in some speci…c situations. In this dynamic context, patients optimally determine a "reservation fee". They will accept all o¤ers charging a fee lower than this reservation fee and keep searching at least one more period if the fee exceeds the reservation fee.

As a consequence of di¤erences in preferences, perceptions, health insurance coverage or income, patients di¤er in their willingness to pay for a medical service (Nay et al., 2016). It goes beyond the purpose of the paper to study the determinants of this willingness to pay. We assume that there are type-H patients characterized by a high willingness to pay H and ( 1) type-L patients with a low willingness to pay L . Let p j be the reservation fee of the type-j patient, with

j 2 [L; H]:
Physicians can chose between practising in the "regulated" sector 1 and in the "free" sector 2.

We assume that the medical service has the same quality regardless of the practitioner. This is a plausible assumption if one takes into account the pure medical service, as a patient should obtain the right treatment regardless of the sector. However, the paper abstracts from the possibility of horizontal di¤erentiation, where physicians in sector 2 can provide better nonessential services, such as better comfort during the waiting time, a reduced waiting time, etc. 8

Let p denote the fee charged by physicians. This fee depends on the sector, the prevailing regulation, and the optimal choice of the physician. In the general framework, we assume that in sector 1, the medical fee is determined by the National Health Insurance in agreement with physicians, and will be denoted p: In France, this regulated fee p is fully reimbursed to the patient by the Social Security and complementary private insurance.

In sector 2, physicians freely determine the fee 9 . Depending on the "target" patients (low or high willingness to pay), they should set as an optimal price either the reservation price p L ; or the reservation price p H : They have no reason to set a fee lower than the reservation fee because patients of a given type would accept the higher reservation fee. The public and private health 8 Such second-order di¤erences can also justify di¤erences in medical fees, in a horizontal di¤erentiation model.

insurances will cover a fraction ( 1) of the di¤erence between the sector 1 fee and the actual fee. In general, > 0; thus patients must pay from their pockets the di¤erence (p p).

The instantaneous utility of a type-j patient who accepts the physician's o¤er p is:

u j (p) = j (p p) ; with j 2 [L; H]; and p 2 [ p; p L ; p H ] (1)
In sector 1, the o¤er is p = p; thus u j ( p) = j :

The production cost of the medical service is normalized to zero. Physicians who chose to work in Sector 2 must bear some additional costs, both …xed (related to degrees and certi…cations) and variable (such as the higher social contributions and the higher representation costs). To keep the model simple, we discard the additional …xed cost and consider that there is only a constant cost t per patient.

Let s denote the frequency of physicians who choose to work in sector 1. Thus, (1 s) is the frequency of physicians who chose to work in sector 2.

Let l denote the frequency of physicians in sector 2 who post the low fee p L ; thus; (1 l) is the frequency of physicians in sector 2 who post the fee p H .

Reservation fees

At each period, patients can …nd either a physician in sector 1 o¤ering the regulated fee p (with frequency s); a physician in sector 2 posting the low fee p L (with the frequency l(1 s)); or a physician in sector 2 posting the high fee p H (with the frequency (1 s)(1 l)):

To determine the reservation fee of the type-L patient, we acknowledge that the patient must be indi¤erent between accepting that o¤er and obtaining the associated instantaneous utility, or rejecting that o¤er and "obtaining" the discounted expected utility from an open sequence of similar searches. Denoting by V L (p) the value function of the type-L patient, his reservation fee (p L ) must verify the indi¤erence condition:

L p L p = sV L ( p) + (1 s) lV L (p L ) + (1 l)V L (p H ) = n s L + (1 s)[l( L p L p ) + (1 l)( L p L p )] o = s 1 (1 s) L : (2)
This result is obtained by acknowledging that the type-L patient who gets an o¤er (p H ) continues to search insofar as the posted fee exceeds his reservation fee. Given the de…nition of the reservation fee, his discounted expected utility from the search activity is equal to ( L p L p ): The implicit de…nition of the reservation fee is:

(p L p) = 1 1 (1 s) L (3)
Obviously, p L > p; the reservation fee of L-type patients is higher than the regulated fee in sector 1.

Following a similar reasoning, and denoting by V H (p) the value function of a type-H patient,

p H p = sV H ( p) + (1 s) lV H (p L ) + (1 l)V H (p H ) = n s H + (1 s) h l( H p L p ) + (1 l)( H p H p ) io = [s + (1 s)l] H l(1 s) p L p 1 (1 s)(1 l) (4)
Replacing p L p by its expression in Eq (3), we obtain:

H p H p = [s + (1 s)l] H l(1 s) 1 1 (1 s) L 1 (1 s)(1 l) or p H p = H [s + (1 s)l] H l(1 s) 1 1 (1 s) L 1 (1 s)(1 l) = 1 1 (1 s)(1 l) H + l(1 s) 1 (1 s) L (5)
Both spreads p L p and p H p are increasing functions in the private coverage rate (1 ):

In other words, if the coverage rate increases, this pushes up the two reservation prices. Because this variable is not essential for the model, in the following we will assume that = 1; i.e., the patient covers from his pocket the whole di¤erence between the sector 2 fee and the sector 1 regulated fee. Then, the two reservation prices have the expressions:

p L = p + 1 1 (1 s) L (6) p H = p + 1 1 (1 s)(1 l) H + l(1 s) 1 (1 s) L (7)

Equilibria

In this game, physicians present three possible pure strategies. We list below these strategies and the associated net payo¤ (y): We recall that t represents the (constant) marginal cost that physicians who choose sector 2 must bear.

Strategy 1 : Work in sector 1 (and charge the regulated fee p). The expected net income associated with this strategy is:

y 1 ( p) = p (8)
Strategy 2 : Work in sector 2 and charge the low fee p L : This o¤er can be accepted by all patients.

The net income is:

y 2 (p L ) = ( p t) + 1 1 (1 s) L : (9)
Strategy 3 : Work in sector 2 and charge the high fee p H : This o¤er can be accepted only by the H type patients, which represent only a proportion of all patients. The expected net income is:

y 3 (p H ) = p H t = ( p t) + (1 ) 1 (1 s)(1 l) H + (1 s)l 1 (1 s) L : (10)
We can show that the game presents three pure strategy equilibria.

Table 1 presents the equilibrium net income and the net income if the physician deviates from the equilibrium strategy, as well as the equilibrium conditions and their range of existence, given

y 1 ( p) p p p y 2 (p L ) y 2 (p L ) = ( p t) + (1 ) L p t + L p t + L y 3 (p H ) h ( p t) + (1 ) H i h ( p t) + (1 ) H + L i p t + H condition 1 (c1) y 1 ( p) > y 2 (p L ) y 2 (p L ) > y 1 ( p) y 3 (p H ) > y 1 ( p) condition 2 (c2) y 1 ( p) > y 3 (p H ) y 2 (p L ) > y 3 (p H ) y 3 (p H ) > y 2 (p L ) c1 ! cost t t> (1 ) L t < L t < H ; otherwise p 3 < 0 c2 ! fee p p > p 1 = 1 h (1 ) H t i p > p 2 = t + [(1 ) H + L ] L (1 ) 8 > > < > > : p < p 4 = (1 ) H t p < p 3 = t + H L (1 ) Equilibrium price p p L = p + L p H = p + H Patient surplus L + ( H L ) H L 0 Table 1.Pure Strategy Equilibria
The last row indicates the patient surplus for each equilibrium, determined as:

W ( p) = (1 ) L + H = L + ( H L ); W (p L ) = (1 )( L p L p ) + H p L p = ( H L ); W (p H ) = (1 )( L p L p ) + H p H p = 0:
If the cost t is too large, either strategy 2 or both strategies 2 and 3 cannot be equilibrium strategies. If t is too small, strategy 1 cannot be an equilibrium strategy. For the sake of parsimony we will analyze only the most general case where all strategies are feasible, which requires the cost condition:

(1

) L < t < L (11)
Table 2 presents the order of critical fees depending on the parameter ; which, we recall, represents the proportion of type-H patients. Condition p 3 < p 4 , t < L (assumed to be ful…lled)

p 2 < p 3 always ful…lled p 1 < p 2 , > 2 = L t L p 1 < p 3 , > 3 = L t H p 1 < p 4 , always ful…lled
Table 2. Order of Critical Fees Obviously, depending on ; we have three cases:

a) p 1 < p 2 , 2 < b) p 2 < p 1 < p 3 , 3 < < 2 c) p 3 < p 1 , < 3 :

Discussion

We now analyze the most plausible case where the number of patients with a high willingness to pay is not very large < 3 = L t H : The other cases are not very di¤erent, they just extend the range of p for which the sector 1 equilibrium is feasible, and thus extends the scope for multiple equilibria. 1. If p < p 2 ; only the equilibrium where all physicians work in sector 2 and post the high fee (p H ) can emerge.

2. If p 2 < p < p 3 ; the game presents a typical con…guration of multiple equilibria; both pure strategy equilibria in which all physicians choose sector 2 (and either set the fee p H or the fee p L ) can occur.

In fact, for p 2 < p < p 3 the game also features a mixed strategy equilibrium where all physicians chose sector 2 (s = 0) yet l 0 choose p L and (1 l 0 ) choose p H : The equilibrium l 0 can be inferred from the condition y 2 (p L js = 0) = y 3 (p H js = 0): After some calculations, its expression is:

l 0 =
(1 ) ( 1)

" H L ( p t) (1 ) L + p t # (12)
We can observe that dl 0 =d p < 0; and that l 0 = 1 for p = p 2 and l 0 = 0 for p = p 3 : If this equilibrium prevails, the higher the regulated fee p; the greater the proportion of physicians (1 l 0 ) that will post the high price p H : 3. If the Health authority aims to make sure that the "bad" equilibrium p H is ruled out for sure, they should impose a fee (slightly) above p = p 3 = t + 4. Policymakers might aim to make sure that the sector 1 equilibrium, ensuring the largest patient surplus, is also feasible. For so doing the government should push the regulated price up to

p = p 1 = 1 h (1
) H t i : Such a higher regulated fee is not innocuous, as it would increase the total budget of the National Health Insurance, and the taxes needed to …nance it.

5. Unfortunately, for p p 1 both equilibria that serve the whole patient market are feasible and, on pure theoretical grounds, it is di¢ cult to state which of the two will prevail. One intuitively appealing re…nement is Pareto-dominance, where we can assume that physicians will prefer the equilibrium where they obtain the highest pro…t. If all physicians choose sector 1, their pro…t margin is y 1 ( p) = p; if all work in sector 2 and post p L their pro…t margin is y 2 (p L ) = p L t = p + L t: Because L t > 0; this criterion for equilibrium selection suggests that the sector 2 equilibrium is most likely to prevail.

In general, systems with multiple equilibria are highly unstable, as changes in patient beliefs, themselves related to the economic outlook of the sector or to policy reforms, can prompt physicians to move from one strategy to another. Our simple model cannot capture the dynamic of change. Furthermore, the assumption according to which all physicians present the same cost of moving from sector 1 to sector 2 explains why pure strategy equilibria are prevalent. Should we relax this assumption and allow the cost t to follow a standard statistical distribution, the same logic would allow the emergence of equilibria with a proportion 0 < s < 1 of physicians choosing sector 1.

6. Finally, we can use the model to analyze policy reforms such as the OPTAM agreement as presented in the Introduction. Physicians who accept this contract agree to maintain fees to maximum 200% of the sector 1 fees. It is realistic to assume that the regulation sets a cap on the highest fees in the market. More precisely, let us denote the new regulated price in sector 2 by p ; with p L < p < p H : The National Health Insurance agreed that complementary private insurance schemes fully cover this fee, while they can cover only (p H p) if the physician applies p H : We have seen that without the cap, p H will be accepted only by the "premium" patients (with a higher willingness to pay). Because the fee is fully reimbursed, now all patients can accept this o¤er, which makes it extremely interesting for the physician. However, there is a big catch in this measure, because no physician has the incentive to post p L : After the adoption of the measure in 2016, there is a documented trend among specialists in sector 2 to push up their fees, which can be interpreted in light of the model as a move from p L to p :

Conclusion

The market for medical services in France has an original two-sector structure. Physicians working in sector 1 must charge a "regulated" fee determined by centralized negotiations between the National Health Insurance and physicians'unions, while physicians working in sector 2 can freely decide on their fees. Specialized physicians can choose the sector in which they perform their medical practice. Potential patients can freely choose their physicians, and when needed, they undertake a standard search process for the best o¤er. This paper introduces a simple equilibrium search model to analyze the dual market for specialist physicians in France. Search models explicitly take into account the cost of producing price information in decentralized markets. Our analysis revealed that this market organization is prone to multiple equilibria with the inherent instability related to these analytical structures. In this context, changes in patients'beliefs can result in large swings of physicians from one sector to another, and even changes from a low to a high fee equilibrium.

We studied in depth the case where the proportion of patients with a high willingness to pay is relatively small. The analysis shows that there is a regulated price that rules out the most ine¢ cient equilibrium where all physicians charge a very high reservation price. For this particular regulated fee, physicians still make positive pro…ts as they prefer to work in sector 2 (but post a low fee). If the Health authority wants to make possible the best equilibrium from the patients' perspective where all physicians work in sector 1, then it must raise the sector 1 regulated fee even further. Unfortunately, these higher fees do not guarantee the emergence of sector 1 equilibrium.

The steadily increasing proportion of physicians opting for sector 2 as observed over the last thirty years in France can be interpreted as a slow transition from the sector 1 equilibrium to the sector 2 equilibrium, inherently reducing patient welfare. In light of our analysis, raising the regulated fee in sector 1 might not reverse this trend. Furthermore, caps on the highest fees posted in sector 2 with full reimbursement as recently enforced by the French Health authorities can also be an ine¢ cient measure, as it would prompt all physicians to choose this fees, including those who initially favoured the lower fees.

Acknowledgement 1 The authors are grateful to Gérard de Pouvourville for his suggestions and

FigureFigure 1 :

 1 Figure1presents the physicians'decision tree:

Figure 2

 2 Figure2presents the regioning of equilibria.

Figure 2 :

 2 Figure 2: The Equilibrium Regioning (for < 3 )

  market power in a market with trade frictions. The larger the frequency ; the higher the margin. Patient welfare is H L :

In 2004, a reform of complementary insurance contributed to further extend the coverage. See Saliba and Ventelou 2007 for an analysis of the reform.

See http://www.leciss.org/newsletter/test-selection-contenu-6

Can be translated as the "option to post moderate fees".

Traditional search literature has acknowledged that the production of price information in decentralized markets is a time consuming, costly phenomenon. See for instance: Stigler, 1960;[START_REF] Mccall | Economics of information and job search[END_REF][START_REF] Diamond | A model of price adjustment[END_REF][START_REF] Diamond | Consumer di¤erences and prices in a search model[END_REF]

0 For p 3 < p < p 1 only the equilibrium where physicians charge the fee p L = p + L exists.


Dmitry V Sorokin
email: sorokin@fi.muni.cz

Igor Peterlik

Marco Tektonidis

Karl Rohr

Pavel Matula

Non-rigid Contour-Based Registration of Cell Nuclei in 2D Live Cell Microscopy Images Using a Dynamic Elasticity Model

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Abstract-The analysis of the pure motion of subnuclear structures without influence of the cell nucleus motion and deformation is essential in live cell imaging. In this work, we propose a 2D contour-based image registration approach for compensation of nucleus motion and deformation. The proposed approach is the extension of our previously published contour-based approach that used an elasticity model to register cell images. Compared to the previous contour-based approach, our approach employs an dynamic elasticity model for forward simulation of nucleus motion and deformation based on the motion of its contours. The contour matching process is embedded as a constraint into the system of equations describing the elastic behavior of the nucleus. This results better performance both in terms of quality of the results and computation time. Our approach was successfully applied to real live cell microscopy image sequences of different types of cells including the data specifically designed to evaluate cell image registration methods. An experimental comparison with existing contour-based registration methods and an intensity-based registration method has been performed. We also studied the dependence of the results on the choice of method parameters.

Index Terms-biomedical image analysis, microscopy, image sequence analysis, registration, elasticity, contour-based registration

I. INTRODUCTION

T HE analysis of subcellular foci motion is one in the most common tasks of live cell imaging. This problem arise in studies related to DNA repair, nucleoli assembly, viral defense and other processes that are important for understanding cell function. However, the motion of foci cannot be sufficiently analyzed directly in raw image sequence as living cells are moving and deforming during imaging process. The observed motion of subcellular foci consists of two components: local motion of the foci and global motion of the nucleus. To get the information about only the subcellular foci motion, the global motion of the nucleus needs to be compensated. This is usually done by means of image registration. Namely all the images in the sequence are normalized to some reference time point (usually the first image of the sequence).

A number of different registration approaches for cell deformation and motion compensation have been presented in previous work. The approaches that use image intensity have been mostly used. In [START_REF] Wilson | A correlation-based approach to calculate rotation and translation of moving cells[END_REF] a rigid approach based on phase correlation has been presented. Another correlationbased approach to compute image translations was proposed in [START_REF] Goobic | Image stabilization and registration for tracking cells in the microvasculature[END_REF]. Baheerathan et al. [START_REF] Baheerathan | Registration of serial sections of mouse liver cell nuclei[END_REF] used phase correlation combined with landmark based method to compensate for affine transformation. The approach presented in [START_REF] Rieger | Alignment of the cell nucleus from labeled proteins only for 4d in vivo imaging[END_REF] used inertia tensor computed based on labeled proteins to compensate for nuclei rotation and translation. In [START_REF] Van De Giessen | Fully automated attenuation measurement and motion correction in flip image sequences[END_REF] the authors used photobleaching model to compensate rigid motion of live cells during FLIP experiments. Raza et al. [START_REF] Raza | Ramtab: Robust alignment of multi-tag bioimages[END_REF] used a block-matching approach to compensate for translations between color channels in multitag fluorescence microscopy images. Besides rigid methods several non-rigid intensity-based registration approaches were presented. In [START_REF] Ozere | Robust parametric stabilization of moving cells with intensity correction in light microscopy image sequences[END_REF] the authors suggested a parametric motion model for cell motion compensation. In [START_REF] Mattes | Analyzing motion and deformation of the cell nucleus for studying co-localizations of nuclear structures[END_REF] the a semiautomatic landmark-based non-rigid registration approach was proposed. In [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] the authors introduced an intensity-based nonrigid registration approach for compensation of cell nuclei motion and deformation based on the extension of Lucas-Kanade optic flow algorithm. The multi-frame extension of this algorithm was proposed in [START_REF] Tektonidis | Non-rigid multi-frame registration of cell nuclei in live cell fluorescence microscopy image data[END_REF]. In [START_REF] Medioni | Motion compensation in two-photon microscopy temporal series[END_REF] the authors suggested a landmark-based approach for affine registration of moving cells. A few algorithms for registration of 2D slices of 3D static microscopy images were introduced. The algorithms based on using thin-plate splines [START_REF] Du | 3d image stack reconstruction in live cell microscopy of drosophila muscles and its validation[END_REF] and B-splines [START_REF] Hogrebe | Serial section registration of axonal confocal microscopy datasets for longrange neural circuit reconstruction[END_REF], [START_REF] Lorenz | A multi-resolution approach to non-rigid registration of microscopy images[END_REF] were presented. The method proposed in [START_REF] Saalfeld | Elastic volume reconstruction from series of ultra-thin microscopy sections[END_REF] was based on block-matching algorithm. The authors in [START_REF] Akselrod-Ballin | Accelerating image registration with the johnson-lindenstrauss lemma: Application to imaging 3-d neural ultrastructure with electron microscopy[END_REF] presented a feature-based registration approach for 3D reconstruction of electron microscopy images. In order to register static confocal microscopy images of different objects an approach based on B-splines was presented in [START_REF] Tomer | Profiling by image registration reveals common origin of annelid mushroom bodies and vertebrate pallium[END_REF] and a landmark-based algorithm that use thin-plate splines was proposed in [START_REF] Peng | Brainaligner: 3d registration atlases of drosophila brains[END_REF].

Intensity-based registration approaches demonstrate high efficiency due to the fact that they use all the information presented in the image sequence. However these methods are not very suitable for image data where the image intensities are not very reliable to be used for image registration (e.g., high level of noise, large and abrupt motion of cellular structures). For example, the methods based on optic flow principles [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF], [10] or block-matching [START_REF] Saalfeld | Elastic volume reconstruction from series of ultra-thin microscopy sections[END_REF] require image intensities to be stable which can not be fulfilled in case of strong intensity changes (e.g. due to photobleaching and noise). In [START_REF] Ozere | Robust parametric stabilization of moving cells with intensity correction in light microscopy image sequences[END_REF] a solution to intensity variation problem was suggested however the proposed registration approach employed affine deformation model which cannot cope with strong local deformations. The landmark-based algorithms [START_REF] Mattes | Analyzing motion and deformation of the cell nucleus for studying co-localizations of nuclear structures[END_REF], [START_REF] Medioni | Motion compensation in two-photon microscopy temporal series[END_REF]- [START_REF] Hogrebe | Serial section registration of axonal confocal microscopy datasets for longrange neural circuit reconstruction[END_REF], [START_REF] Akselrod-Ballin | Accelerating image registration with the johnson-lindenstrauss lemma: Application to imaging 3-d neural ultrastructure with electron microscopy[END_REF], [START_REF] Peng | Brainaligner: 3d registration atlases of drosophila brains[END_REF] require stable landmarks to be localized in the registered pair of images which can be difficult as often there is no a priori biologically stable structures in the moving cell. Also, when the subject of the study is the local motion of subcellular structures the intensity-based registration algorithms can influence the motion patterns of these structures while compensating the cell motion and deformation. Thus the errors can be introduced into the trajectories of the structures of interest. This problem can be solved by acquiring the images in 2 channels [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF], [START_REF] Tektonidis | Non-rigid multi-frame registration of cell nuclei in live cell fluorescence microscopy image data[END_REF] where the objects of interest are represented in one channel and some stable structure (e.g. chromatin) is represented in the second channel. The approaches [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF], [START_REF] Tektonidis | Non-rigid multi-frame registration of cell nuclei in live cell fluorescence microscopy image data[END_REF] were applied to register the nucleus channel of two-channel microscopy image sequences and the computed transformations could be used to compensate the nucleus motion and deformation in the particle channel. However it is not always possible to acquire multichannel images due to different biological experimental setups. The mentioned limitations of intensitybased registration approaches led to the need of development of cell motion compensation approaches that do not use image intensity for registration.

In previous work on registration of cell images in microscopy data, only a few approaches that do not use image intensity have been presented. In [START_REF] Sage | Automatic tracking of individual fluorescence particles: application to the study of chromosome dynamics[END_REF] the authors used a modelbased approach where yeast cells were fitted with an ellipses to compensate for rigid motion in 2D image sequences. Matula et al. [START_REF] Matula | Fast point-based 3-d alignment of live cells[END_REF] introduced a point-based rigid registration approach where the positions of subcellular structures were used to find the suitable transformation between the cell in different frames. A shape based approach using distance functions was introduced in [START_REF] Dzyubachyk | Automated analysis of time-lapse fluorescence microscopy images: from live cell images to intracellular foci[END_REF]. Several non-rigid registration approaches without using image intensity were proposed. In [START_REF] Yang | Nonrigid registration of 3-d multichannel microscopy images of cell nuclei[END_REF] the authors used an extension of demons algorithm applied to cell segmentation masks to compensate for global motion and deformation of living cells. De Vylder et al. [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] introduced a contour-based non-rigid registration method for dynamic cell nuclei where the deformation is defined based on the correspondences between the contour points of the nucleus in consecutive time frames.

The described approaches do not use image intensity but binary segmentations of the nucleus. The efficiency of the methods were demonstrated on both real microscopy data and synthetically created cell phantom images. These approaches were mostly applied to the data acquired in one channel. However in [START_REF] Sage | Automatic tracking of individual fluorescence particles: application to the study of chromosome dynamics[END_REF], [START_REF] Matula | Fast point-based 3-d alignment of live cells[END_REF] the algorithms could compensate only the rigid motion of the cell. In [START_REF] Dzyubachyk | Automated analysis of time-lapse fluorescence microscopy images: from live cell images to intracellular foci[END_REF], [START_REF] Yang | Nonrigid registration of 3-d multichannel microscopy images of cell nuclei[END_REF] the authors described non-rigid approaches but the implied deformation models allow to compensate only limited deformations that is not very effective in the case of abrupt motion. In [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] the approach was specifically designed to cope with strong deformations however due to the uniform sampling of the nucleus boundary the contour matching algorithm introduces errors in the deformation fields. In addition, the thin-plate deformation model used in [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] is a relatively rough approximation of the nucleus deformation. In [START_REF] Gladilin | Shape normalization of 3d cell nuclei using elastic spherical mapping[END_REF] the authors introduced a cell nucleus shape normalization approach based on linear elasticity. However the approach [START_REF] Gladilin | Shape normalization of 3d cell nuclei using elastic spherical mapping[END_REF] was applied to static cell nuclei to normalize their shape to 3D sphere and measure relative position of internal structures.

In this contribution, we propose a new contour-based nonrigid registration approach. In our approach we use dynamic non-linear elasticity model to compensate for cell motion and deformation. This work extends our previous work [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] where elastostatics model employing Navier equation with Dirichlet boundary conditions was used to register the cell nuclei. The method described in [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] employed "global-tolocal" registration strategy that required two steps to be performed: the compensation for nucleus motion by rigid alignment of the images and the compensation for nucleus deformation by solving the Navier equation. The registration was performed for each pair of two consecutive frames resulting the backward deformation field from frame k to frame k -1. The Navier equation was solved using the Dirichlet boundary conditions defined by the contour matching of two consecutive frames. As distinct from the previous approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF], the method presented in this contribution is based on the forward simulation of the elastic motion of the nucleus resulting the full model of nucleus deformation given on triangulation that can be used to construct both forward and backward deformation fields for cell image registration. The proposed approach does not require rigid alignment step as the new contour matching procedure is capable to compensate for large rigid displacements. The contour matching is embedded as constraints into the system of equations of the elastic motion model. The system of equations is iteratively solved using prediction-correction scheme. We also significantly extended the quantitative evaluation of the proposed approach. The method was evaluated on real fluorescent microscopy image data including the data specifically designed and acquired for testing cell image registration techniques [START_REF] Sorokin | Visualizing stable features in live cell nucleus for evaluation of the cell global motion compensation[END_REF]. The approach was compared with the method described in our previous work [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF], the contour-based method presented in [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] and the intensity-based method [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF].

The paper is organized as follows. In the next section we introduce our contour-based non-rigid registration approach. Then we present the evaluation results and method parameter study. In the final sections of the paper we conclude the paper with a discussion.

II. METHODS

In this section, we describe our contour-based registration method. Prior to registration we perform threshold-based segmentation and determine the contours of cell nuclei. Below, we first describe the data representation and overall scheme of our registration approach (see Fig. 2). Then, we explain the physical model used for the registration and its difference from our previous approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF]. After that we detail the main parts of the approach, namely the simulation of nucleus motion, the definition of contour matching constraints, and the numerical solution scheme.

A. Data Representation and Forward Nucleus Motion Simulation

To perform the registration of every image of the sequence to the first image we use the forward simulation of nucleus motion based on the model described in the next subsection based on finite elements (FE) method. Given the binary images of cell nucleus in every frame we define nucleus contours C k where k = 0, ..., M -1 and M is the number of images in the sequence. Each point of the contour is represented by its position P and outer normal n P . Then we define the triangulation T 0 with N nodes for the first frame of the sequence where E 0 is the set of triangulation elements and N 0 is the set of nodes. The idea of forward simulation of nucleus motion is to simulate the elastic motion of the nucleus given only the contours C k . In order to do this, the triangulation T 0 is iteratively deformed based on the model described below so that its boundary optimally fits the contour C 1 . Let T 1 be the result of this operation representing the deformed nucleus in the second frame of the sequence. The same operation is performed consequently for T k and C k+1 resulting the set of triangulations T k (k = 1, ..., M -1) that represent the deformed nucleus in the corresponding frame.

Given the triangulation representing the deformed nucleus in every frame T k (k = 0, ..., M -1) one can compute the forward deformation field u 0→k = q kq 0 or backward deformation field u 0←k = q 0 -q k where q 0 , q k ∈ R 2N are the vectors that consist of T 0 and T k triangulations node positions (q = [x 0 , y 0 , x 1 , y 1 , ..., x N -1 , y N -1 ] T ). Thus u 0→k and u 0←k are defined in the nodes of triangulation. The forward deformation fields u 0→k can be interpolated in every pixel of the image and used to wrap the image sequence to normalize it to the first frame of the sequence. The backward deformation field u 0←k can be used to register a set of features (e.g. the position of foci the local motion of which is studied).

B. Modeling the Nucleus Deformation Based on Linearized Elasticity and Dynamic Integration

As it was shown in [START_REF] Tseng | Microorganization and visco-elasticity of the interphase nucleus revealed by particle nanotracking[END_REF] cell nucleus demonstrates elastic behavior. Thus, the deformation of the nucleus can be modeled using elasticity theory [START_REF] Sadd | Elasticity: Theory, Applications, and Numerics[END_REF] considering cell nucleus an isotropic homogenous material.

The nucleus motion can be described using the second Newton's law

M a = f int + f , f int = µ∆u + (µ + λ)∇(∇ • u), (1)
where M is the nucleus mass, a is the acceleration, u represents the displacement and f represents the external forces applied to the nucleus boundary due to constrains that drive the registration. The Lamé coefficients µ and λ define the properties of the elastic material and connected to the Poisson's ratio and Young's modulus with the following relations

λ = νE (1 + ν)(1 -2ν) , µ = E 2(1 + ν) .
The f int term represents the internal elastic forces that present in the nucleus. In [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] we used the elastostatic approach exploiting Navier-Cauchy equation to model the nucleus motion:

f int = 0,
where f int is defined in (1) and the contour matching constrains that drove the deformation were defined as Dirichlet boundary conditions. In this contribution we employ the dynamic model by adding the M a term responsible for the nucleus inertia. The contour matching constrains that drive the deformation are defined using the external forces f . In order to solve (1) we employ FE method. Let us consider the nucleus triangulation T k in frame k and the target contour of the nucleus C k+1 in frame k + 1. Then the motion of the nucleus can be presented by the Newton's second law rewritten in the following form

Mq = F(q, q), H(q, p) = 0, (2)
where q ∈ R 2N is the vector of triangulation node positions, p ∈ R 2Nc is the vector of points in the target contour C k+1 (N c is the number of points in C k+1 ), M ∈ R 2N ×2N is the mass matrix (the nucleus is considered uniform and thus the mass matrix does not depend on q). The internal elastic forces component f int in ( 1) is represented with the mapping F(q, q) : R 2N → R 2N . The constraints given by the external forces f in (1) on the boundary edges of T k are given in the implicit form using the mapping H(q, p) : R 2Nc → R 2Nc . The problem ( 2) is considered as a constrained optimization problem. We use the Lagrange multipliers method to solve it

Mq = F(q, q) + H T (q)λ, (3)
where H(q) = ∂H ∂q is the Jacobian of the constraint mapping H(q, p) and λ ∈ R 2Nc is the corresponding vector of Lagrange multipliers. Thus H T (q) defines the direction of constraint forces and λ represents the constraint force amplitudes (the details of constraints definition is described in section II-C).

Let us consider the time interval [t k , t k + h]. Integrating (3) with respect to time t we get

M( qk+h -qk ) = t k +h t k F(q, q) dt + hH T (q)λ, q k+h = q k + t k +h t k q dt. (4)
We use the implicit Euler method to compute the integrals in the right part of Eq. ( 4): M( qk+h -qk )=hF(q k+h , qk+h ) + hH T (q k+h )λ k+h ,

q k+h =q k + h qk+h . ( 5
)
As it is described in section II-B we assume that the deformations are small thus we can linearize F using Taylor expansion:

F(q k + dq, qk + d q) = f k + δF δq dq + δF δ q d q. ( 6
)
Let us denote K = δF δq and B = δF δ q . These matrices are known as the stiffness and damping matrices. The details on K and B matrices construction can be found in [START_REF] Nesme | Efficient, Physically Plausible Finite Elements[END_REF]. In order to compensate for effects appearing in the case of large displacements the corrotational formulation for stiffness matrix construction is used [START_REF] Felippa | A study of optimal membrane triangles with drilling freedoms[END_REF]. The explicit form of these matrices is given in Appendix A. Substituting (6) into (5) and using dq = q k+h -q k = h qk+h and d q = qk+h -qk we obtain

(M + hB + h 2 K) A d q = -h 2 K qk b +hH T (q k+h )λ k+h . (7)
This equation defines the change in velocity in the nodes of the nucleus triangulation during one time step. The method parameters that control the elastic properties of the nucleus, namely Young's modulus E and Poisson's ratio ν, are contained in the matrix K. The only unknown term in this equation is H(q k+h )λ k+h the computation of which is discussed in the next sections.

C. Definition of Contour Matching Constraints

Let us consider the nucleus triangulation T k in frame k and target contour of the nucleus C k+1 in frame k + 1. Let us also denote the subset of edges T k that connects the boundary points of the triangulation T k (see Fig. 3a). For every point P ∈ C k+1 we find the closest point Q c ∈ T k that n P • n Q > 0 where n P and n Q are the normals in points P and Q c (Q c and n Q are calculated as linear interpolation of neighboring nodes of T k ). The latter condition is necessary to avoid points Q c located on the opposite side of the nucleus in the case of large frame displacements between frames (see Fig. 3b). Then we define the gap between P and Q c as

δ P = (P -Q c ) • n Q . ( 8
)
This gap is zero when the point P is located on the tangent of T k in the point Q c . Such gaps are defined for all P i ∈ T k , i = 0, N c -1. The corresponding points Q ci ∈ T k are obtained using the ICP algorithm [START_REF] Rusinkiewicz | Efficient variants of the icp algorithm[END_REF]. We consider the contours to be matched when δ Pi = 0 for ∀i = 0, N c -1. Thus the constraint mapping H(q, p) can be defined as

H(q, p) =    δ P0
. . .

δ P Nc -1    , P i ∈ C k+1 .
It is important to note that the points Q ci ∈ T k are the points of the same elastic body. Thus applying the constraint in one point influence the other points. So the optimal minimization of the gaps (i.e. contour matching) taking into account the elasticity of the body is obtained using the Lagrange multipliers λ k+h in [START_REF] Ozere | Robust parametric stabilization of moving cells with intensity correction in light microscopy image sequences[END_REF]. To define λ k+h and solve [START_REF] Ozere | Robust parametric stabilization of moving cells with intensity correction in light microscopy image sequences[END_REF] we use a prediction-correction solution scheme which is described in the next subsection.

D. Prediction-Correction Solution Scheme

In this section we describe the computation scheme allowing to solve the equation [START_REF] Ozere | Robust parametric stabilization of moving cells with intensity correction in light microscopy image sequences[END_REF]. First, to simplify the solution process we assume that the constraint Jacobian H T (q) does not change during constraints resolving thus instead of H T (q k+h ) we use H T (q k ) = H T . In the following, we use the notation H to emphasize that this matrix does not change during time step. Then, substituting this into (7) we obtain

A d q = b + hH T λ k+h . (9)
The solution consists of three steps: prediction, constraint solving and correction. Let us describe these steps.

Step 1. First we solve (9) setting λ k+h = 0

A d qp = b. ( 10
)
Thus we obtain the predictive change in velocity d qp that corresponds to the unconstrained motion that continue from the previous step (analogous to body inertia).

Step 2. On this step we resolve the constraints and get the correction for the vector of Lagrange multipliers λ k+h . Assuming b = 0 we get

A d qc = hH T λ k+h . (11)
In [START_REF] Medioni | Motion compensation in two-photon microscopy temporal series[END_REF] the correction change in velocity d qc and the vector of Lagrange multipliers λ k+h are unknown. To obtain the relation between these variables we linearize the constraints equation with respect to velocity change:

H(q k+h , p) = H(q k , p) + hH T d qc . ( 12
)
According to (2) the constraints are satisfied when H(q k+h , p) = 0. Taking this into account and substituting the solution of ( 11) into [START_REF] Du | 3d image stack reconstruction in live cell microscopy of drosophila muscles and its validation[END_REF] we get the following non-linear equation

H(q k , p) + hHA -1 H T λ k+h = 0. (13)
In this equation the term C q = HA -1 H T called source object compliance matrix is responsible for the elastic behavior of the nucleus in the constraints solution. In [START_REF] Hogrebe | Serial section registration of axonal confocal microscopy datasets for longrange neural circuit reconstruction[END_REF] any physical properties of the target object represented by C k+1 are not taken into account as this object is a purely geometrical entity in this simulation. Although the solution of (13) can suffer from inaccuracies in extraction of C k+1 from the source images. To increase the numerical stability of the solution process we introduce to (13) the target object compliance matrix C p similar to [START_REF] Courtecuisse | Real-time simulation of contact and cutting of heterogeneous soft-tissues[END_REF], [START_REF] Duriez | Realistic haptic rendering of interacting deformable objects in virtual environments[END_REF]. However since the object C k+1 does not have any physical properties, the matrix C p is diagonal: C p = γI 2Nc . The additional method parameter γ is called target compliance. Thus we obtain the following equation

H(q k , p) + h(C q + C p )λ k+h = 0. ( 14
)
which is solved by Gauss-Seidel method with respect to λ k+h .

Step 3. After we obtain the values of λ k+h we can compute the correction change in velocity d qc using [START_REF] Medioni | Motion compensation in two-photon microscopy temporal series[END_REF]. We compute the overall change in velocity during time step as

d q = d qp + d qc = A -1 b + h 2 A -1 H T λ k+h . (15)
Then we can obtain the velocity on the next time step as qk+h = qk + d q and subsequently the triangulation nodes positions on the next time step q k+h = q k + h qk+h . Thus the problem ( 5) is solved and we can proceed to the next time step.

The iterative solution process is continued until the equilibrium is reached i.e. the triangulation T k representing the nucleus in frame k fits the contour C k+1 . Typically 20 time steps is enough to reach the equilibrium. Then we proceed to the next pair of frames taking the deformed triangulation as the input. The overall scheme of the simulation process is shown in Fig. 2.

III. EXPERIMENTAL RESULTS

For quantitative evaluation of our approach we used two datasets of real microscopy image sequences of live cells. The image sequences in both datasets contained the structures that are stable w.r.t. the nucleus (i.e., there is no relative motion). These stable structures were used to assess the registration error of the proposed algorithm in comparison to the contourbased approaches from [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF], [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] and intensity-based approach from [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF]. The images from the first dataset contain local stable features manually tracked by human annotators that allowed to study the methods in case of significant rotation, translation and shrinkage of the cell. The images from the second dataset contain automatically tracked features that occupy significant part of the nucleus and enable to compare the methods efficiency in the interior part of the nucleus and on its periphery. The sequences from the second dataset also have strong intensity variations and nucleus deformation which allows to test the methods in the case of intensity changes.

A. Evaluation Datasets Description

The first dataset consists of three image sequences that were used for the evaluation of the approach in [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF]. In these image sequences the stable structures were manually tracked by annotators. In two sequences (denoted as Seq. A1 and Seq. A2) the nucleus of U2OS cell was UV-irradiated in a stripe-like region. The UV-irradiated stripe is visible in the nucleus channel and stable w.r.t. the nucleus. In both sequences keypoints of the UV-irradiated stripe were manually tracked by three independent annotators. The nuclei in these image sequences undergo significant translation (including jumps) and rotation along with nucleus deformations. In addition, the image noise is relatively high.

The third image sequence (denoted as Seq. A3) depicts U2OS cell nucleus going into mitosis. For the evaluation 9 spot-like structures in the nucleus channel are used as keypoints that were manually tracked by one annotator. The structures are stable w.r.t. the nucleus and represent the motion and deformation of the nucleus. For Seq. A1 and Seq. A2 the ground truth was obtained by manual tracking by three different annotators to increase the reliability of the evaluation, since the data is quite noisy and the structures are rather blured. In Seq. A3 the spot-like structures are better visible and thus could be tracked reliably by one annotator.

The second dataset contains four image sequences (denoted as Seq. B1-B4) obtained using procedure described in [START_REF] Sorokin | Visualizing stable features in live cell nucleus for evaluation of the cell global motion compensation[END_REF]. These sequences depict HeLa cells with histone H2B tagged by GFP. In each sequence four regions of interest were bleached using argon laser. The bleached structures represent a grid of four orthogonally crossing lines (see Fig. 4). These bleached structures are stable w.r.t. the nucleus and represent the motion and deformation of the nucleus [START_REF] Sorokin | Visualizing stable features in live cell nucleus for evaluation of the cell global motion compensation[END_REF]. In all the sequences the bleached structures were automatically detected using the algorithm described in [START_REF] Kondrat | Automatic detection of laserinduced structures in live cell fluorescent microscopy images using snakes with geometric constraints[END_REF] resulting line features that were used for the evaluation of the registration approach. In addition, we studied the behavior of the registration approaches on the periphery of the nucleus using the end points of the line features (further referred as end points), and on the interior part of the nucleus using the intersections of the line features (further referred as inner points, see Fig. 4).

The details of the image sequences A1-A3, B1-B4 are given in Table I.

B. Evaluation Results

For each keypoint in Seq. A1-A3 we computed the registration error as the Euclidean distance to the position of the keypoint in the first image: where p k i is the coordinate of the i-th feature point in the k-th frame. The same registration error was computed for the inner points and end points in Seq. B1-B4.

e k i = p k i -p 1 i ,
For line-features in Seq. B1-B4 we computed the registration error for every frame as the Fréchet distance to the line-features in the first frame. The parameters for the tested approaches were set to the optimal values to suit our data. For the contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] we used 100 contour sampling points and thin-plate splines for the deformation field interpolation. For contour matching we used the normalized centroid distance signature in Seq. A1, A2, B2, the curvature signature in Seq. A3, B1 and the tangent direction signature in Seq. B3, B4. For the intensity-based approach [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] we used symmetric weighting scheme and σ = 2 for regularization of the deformation and update field. For both elasticity-based approaches we used the parameter values according to the real properties of cell nucleus. The Poisson's ratio of the nucleus was set to 0.4 [START_REF] Gladilin | Contactless determination of nuclear compressibility using 3d image-and model-based analysis of drug-induced cellular deformation[END_REF]. The Young's modulus of the nucleus was set to 10 kPa [START_REF] Lammerding | Mechanics of the nucleus[END_REF] for Seq. A1-A2 and Seq. B1-B4. For Seq. A3 we set the Young's modulus to 2 kPa as the Young's modulus of the nucleus in mitosis was reported to drop 5 times [START_REF] Houchmandzadeh | Elasticity and structure of eukaryote chromosomes studied by micromanipulation and micropipette aspiration[END_REF]. The triangulation for the FE method in both elasticitybased approaches was constructed based on nucleus contours represented by 100 sampling points (to be consistent with the number of contour points used in contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF]). The target compliance for the proposed approach was set to γ = 0.5.

The registration error ēmean averaged over all time points and all keypoints is shown in Table II for Seq. A1 and Seq. A2 for the three different annotators. In addition, the average values over the three annotators has been included. It can be seen that the variation of ēmean for the different annotators for the unregistered case is rather low which indicates the consistency of the selected keypoints. The average values for ēmean over the three annotators show that the proposed approach outperforms the previous contour-based approaches [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] and [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] by 53% and 3% respectively for Seq. A1, and by 127% and 3% respectively for Seq. A2. Compared to the intensity-based approach [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] the proposed approach yields a 27% improvement for Seq. A1 and 10% improvement for Seq. A2 while the intensity-based approach exploits more information using the image intensities. The registration error ēmean for Seq. A3 is also shown in Table II. For Seq. A3 the result of the proposed approach is better compared to the previously introduced contour-based approaches by 48% and 10% respectively but worse compared to the intensity-based approach. The reason is that the intensity-based approach exploits the intensity information within the nucleus (which is going into mitosis) while the contour-based approaches do not depend on the intensity variation which can be an advantage for certain data.

The registration error ēmean averaged over all time points and all line features for Seq. B1-B4 is shown in Table III. The proposed approach gives better results than both previously introduced contour-based approaches for all the sequences from the second dataset. Namely, the proposed approach yields in average 76% improvement in comparison to the contourbased approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] and 13% improvement in comparison to the contour-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF]. In addition, the proposed approach outperforms (in average by 82%) the intensitybased approach for all sequences in the second dataset. The reason for low performance of the intensity-based method for the second dataset is that there is a substantial change in nucleus intensity throughout the sequences caused by imaging conditions. While the method [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] is based on the extension of Lucas-Kanade optic flow algorithm which requires the image intensity to be stable.

We also analyzed the behavior of the proposed approach in the central part of the nucleus and its periphery by using the inner points and the end points of the line features. The comparison of the registration error averaged over all time points for end points and inner points is given in Table III. For inner points the proposed approach yields in average 17% improvement in comparison to the contour-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF], 58% in comparison to the contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] and 103% improvement in comparison to the intensity-based method [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF]. The most significant improvement of the proposed approach in comparison to the previously introduced algorithms was shown for the periphery of the nucleus. The registration error for the end points for the proposed approach was 41% less than for the contour-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF], 104% less than for the intensity based approach [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] and 185% less than the contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF]. However in Seq. B4 the proposed approach performs 17% worse for end points compared to the contour-based method [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF].

To demonstrate the dependence of the registration error on time we show the registration error averaged over keypoints tracked by the annotator A 2 for Seq. A2 (Fig. 6). It can be seen that both elasticity-based approaches have similar behavior and yield better results than the other methods almost through the whole sequence. However the proposed approach performs slightly better in the beginning (time points 1-6) and in the end of the sequence (starting from time point 23). It can also be seen that the error fluctuations of the registration error for the contour-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF], intensity-based approach [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] and the proposed approach is significantly less than the fluctuations for the contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF]. We also show the registration error averaged over line features for Seq. B3 in Fig. 7. Here, the registration error of the proposed approach is the lowest almost through the whole sequence and the behavior of the registration error for all contour-based approaches is similar to Seq. A2. The registration error of the intensity-based approach is higher than in Seq. A2 due to intensity changes. However the error for the proposed approach and the contourbased approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] slightly increases starting from time point 25 while the increase rate for the intensity-based approach is much less.

To visualize the registration error we showed the tracks of the line features overlaid over the first frame of the Seq. B2 for unregistered data and after registration with different approaches. For better visual perception the line features are represented by points sampled with 30 pixel interval. It can be observed that the tracks of the points after registration with both elasticity based approaches are much more compact than after the registration with the contour-based approach [START_REF] Vylder | 2d mapping of strongly deformable cell nuclei-based on contour matching[END_REF] and intensity-based approach [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF] indicating the lower registration error. Moreover, the tracks obtained after the registration using the proposed approach are slightly more compact than after registration using the previous elasticity-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] (the particular tracks are marked with white arrows).

C. Method Parameters Study

We have also studied the dependence of the registration results on the choice of method parameters. As it is described in section II the proposed approach have three parameters: the Poisson's ratio ν, the Young's modulus E and the target compliance γ. We used Seq. A2, A3 and B2 to study the influence of the parameters on the registration error for First, we examined the influence of the Poisson's ratio on the registration results. We performed the computations of ēmean for the fixed values of Young's modulus and target compliance (E = 10kP a, γ = 0.5) and the value of ν ∈ [0.38, 0.49] with the step size 0.01. The range of the tested values of ν was selected according to the latest studies on the Poisson's ratio value for cell nucleus [START_REF] Gladilin | Contactless determination of nuclear compressibility using 3d image-and model-based analysis of drug-induced cellular deformation[END_REF], [START_REF] Avalos | Probing compressibility of the nuclear interior in wild-type and lamin deficient cells using microscopic imaging and computational modeling[END_REF]. For Seq. A2 we took the registration error ēmean averaged over annotators, for Seq. B2 we computed the registration error ēmean for line features. The standard deviation of the registration error ēmean for Seq. A2, Seq. A3 and Seq. B2 were 0.008, 0.011 and 0.007 pixels respectively which is less than 0.2% in relation to the absolute values of the registration error. Thus the value of Poisson's ratio does not significantly influence the registration results if chosen from the range reported in previous work.

Next, we studied the behavior of the registration error depending on the values of E and γ for the fixed value of ν = 0.4. We computed the registration error ēmean of the proposed approach for the fixed value of ν and every pair of E and γ which values were chosen from the selected ranges. Recent studies show that the Young's modulus of the cell nucleus vary from 1 to 10 kPa [START_REF] Lammerding | Mechanics of the nucleus[END_REF], [START_REF] Caille | Contribution of the nucleus to the mechanical properties of endothelial cells[END_REF]. However we extended the range of E for our evaluations and took the values from the interval 0,5 to 15 kPa with the step of 1 kPa. The target compliance regulates the stiffness of the contour matching process by enabling the purely geometrical target object to move to improve the numerical stability of dynamic integration. Namely, the target compliance allows the contour points of the target nucleus to slide along its tangent during the simulation process. Since γ does not have a physical meaning for our evaluation we took γ ∈ [0, 6] with the step of 0.5. The standard deviation of averaged over annotators registration error ēmean for Seq. A2 was 0.033 pixels. The standard deviation of the registration error ēmean for Seq. A3 was 0.011 pixels. And the standard deviation of the registration error computed for line features for Seq. B2 was 0.083 pixels. In Fig. 8 the mean registration error depending on E and γ and the histograms of mean registration error are shown. One can see that for Seq. A3 and Seq. B2 the registration error does not significantly depend on γ but depends on E. However the theoretical values of E and γ (plotted as green dot) taken for the evaluation results presented in Table II and III are within 10 optimal E-γ pairs resulting the lowest registration errors (plotted as red crosses). For Seq. A2 one can see the inverse relationship between the optimal values of E and γ.

The optimal E-γ pairs resulting 10 lowest registration errors are within E ∈ [START_REF] Medioni | Motion compensation in two-photon microscopy temporal series[END_REF][START_REF] Saalfeld | Elastic volume reconstruction from series of ultra-thin microscopy sections[END_REF] kPa and γ ∈ [2, 3] which is higher than the theoretical range 1-10 kPa [START_REF] Lammerding | Mechanics of the nucleus[END_REF]. However the difference between the registration error obtained for E = 10 kPa and γ = 0.5 from the minimal registration error obtained for E = 14 kPa and γ = 2 is just 0.06 pixels.

IV. DISCUSSION

We performed the evaluation of the proposed registration approach on two real microscopy image datasets specifically designed for cell image registration evaluation. It turned out that the proposed approach outperforms previous contourbased methods. The reason for that is the realistic dynamic elasticity model used to simulate the nucleus motion based on the motion of its contours. While in the previously introduced approach [23] any physical model was not employed and in our previously introduced contour-based approach [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] the less accurate elastostatic model was used.

The proposed approach also outperforms pairwise intensitybased registration approach for 6 out of 7 testing sequences. Generally speaking, intensity-based approaches are more beneficial than contour-based approaches as they use the information from entire image while contour-based approaches use only the information about the motion of cell contours. However intensity-based approaches usually does not imply any physical model of motion and rely only on the information taken directly from the images. This can cause problems if the images have high level of noise or significant structural changes in intensity during cell motion. Another limitation of intensity-based approaches is that they require some stable structures to be present in the images to be registered. These structures should have very limited local motion and mostly represent the motion and deformation of the cell. In order to fulfill this requirement the images are usually acquired in two channels where one channel represents the nucleus and another channel represents the particles of interest. Then the intensity-based registration approach is applied to the nucleus channel and the resulting deformation fields are used to compensate the motion of the nucleus in the particle channel [START_REF] Kim | Nonrigid registration of 2-d and 3-d dynamic cell nuclei images for improved classification of subcellular particle motion[END_REF]. However in some cases the image sequence can not be acquired in two channels and both nucleus and particles are shown in the same channel. In this case the application of intensity-based approach can influence the local motion of the particles which is usually the main biological interest in this type of applications. Taking this into account, the contourbased approaches can be beneficial and more universal for data acquired in one channel even though in some cases the registration accuracy of contour-based methods are lower than the accuracy of the intensity-based methods.

In was shown that the proposed approach is quite robust to the choice of method parameter. Two out of three method parameters (Poisson's ratio and Young's modulus) have clear physical meaning and can be chosen based on the theoretical values reported in the literature. Even though there is still no consensus on the exact values of these parameters for different cell phenotypes we have showed that variation of E and ν in the range reported in literature (E ∈ [1, 10] kPa and ν ∈ [0.38, 0.48]) do not lead to significant difference in registration results. The third parameter of the method (target object complience) also does not significantly influence the registration quality and can be chosen empirically.

Unlike the approach in [START_REF] Sorokin | Non-rigid contour-based temporal registration of 2D cell nuclei images using the Navier equation[END_REF] the proposed approach does not require separate rigid alignment step. Thanks to the model where the contour matching is embedded as constraints into the system of equation the proposed approach is capable to compensate for large displacements and rotation of the cell. In addition the solution is based on FE method and consists of three main steps which require solving the linear system of equations, setting the constraints by iterative closest point method and resolving the constraints by Gauss-Seidel method. It makes the calculations computationally efficient. We used SOFA framework to compute the simulation [START_REF] Allard | Sofa-an open source framework for medical simulation[END_REF] and Matlab to compute the deformation fields. We determined the computation time for Seq. A2 containing 38 frames of 512x512 pixels on Windows workstation with an Intel Core 2 Quad CPU (2.4 GHz). The computation time was 3 min 46 sec. Note, that our implementation is not optimized and the computation time can be further reduced.

V. CONCLUSION

We have presented a contour-based approach for temporal non-rigid registration of cell nuclei in 2D live cell microscopy images. A quantitative experimental evaluation of the approach using real live cell microscopy data has been performed. We also investigated the influence of method parameters on registration results. The presented approach exploits forward nucleus motion simulation computed based on FE method. The resulting motion model defined on the set of triangulations allows computing the deformation fields to perform the image registration. By using dynamic elasticity model with contour matching embedded into the system of equation as constraints the proposed approach outperforms the existing contour-based approaches for all tested image sequences and pairwise intensity-based approach for all except one tested image sequences. Thanks to FE method used for calculations the proposed approach is computationally effective. In comparison to the existing contour-based approaches our approach can be straightforwardly extended to 3D by choosing the appropriate finite elements which is a topic for future work. In addition, the proposed solution scheme allows embedding some intensity information to the elastic deformation model by including additional constraints to the equations which is also planned for future work.

APPENDIX A

The stiffness matrix K is constructed by summing up the local stiffness matrices K e for for every element e ∈ E k of the triangulation T k (this process is called the assembly [START_REF] Baraff | Large steps in cloth simulation[END_REF]). The local stiffness matrix for e element is defined as and we refer the reader to [START_REF] Nesme | Efficient, Physically Plausible Finite Elements[END_REF] for the definition of the matrix R e . The damping matrix B is defined as

K e = R
B = αK + βM,
where α and β are the Rayleigh coefficients (in our approach we use α = 10 and β =?), M is the nucleus mass matrix and K is the stiffness matrix.

  (a) Two-channel data. Structures of interest presented in green channel. Nucleus presented in red channel.

Fig. 1 .

 1 Fig. 1. Example images from the time-lapse image sequences (contrast-enhanced).

Fig. 2 .

 2 Fig. 2. Forward nucleus motion modeling scheme.

Fig. 3 .

 3 Fig. 3. (a) Constrain forces definition scheme. The green polygon line represents the nucleus contour T k in frame k. The red polygon line represents the target nucleus contour C k+1 in frame k + 1. The orange arrow represents the constraint force δ P applied to the point Qc. Black dotted lines represent the tangent of T k in point Qc and a parallel to this tangent in the point P . (b). The case when np • nq < 0.

Fig. 4 .

 4 Fig. 4. The first image of the image sequences B1-B4 from the second evaluation detaset. The red lines represent the line features, the green dots represent the inner points, and the yellow dots represent the end points. The images are cropped and scaled to fit the figure.

Fig. 5 .

 5 Fig. 5. Tracks of line features of Seq. B2 overlaid with the first image of the sequence. The tracks of the line features are represented as the tracks of the points of the line features sampled with 30 pixel interval for better visibility. The tracks of line features are shown for (a) unregistered data, (b) after registration performed with the contour-based approach [23], (c) after registration performed with the intensity-based approach [9], (d) after registration performed with the contour-based approach [25], and (e) after registration performed with the proposed approach. White arrows indicate the tracks with the most visible difference between (d) and (e).

Fig. 6 .

 6 Fig. 6. Average registration error for the structures manually tracked by annotator A 1 in the image sequence Seq. A2 as a function of time.

Fig. 7 .

 7 Fig. 7. Average registration error for the line features in the image sequence Seq. B3 as a function of time.

Fig. 8 .

 8 Fig. 8. The mean registration error of the proposed approach depending on E and γ and the histograms of mean registration error for different pairs of E-γ parameters and fixed value of Poisson's ration ν = 0.4. In (a) the mean registration error averaged over three annotators for Seq. A2 is shown, in (b) the mean registration error for Seq. A3 is shown, in (c) the mean registration error for line features for Seq. B2 is shown. The red crosses represent the parameter pairs for which the 10 lowest mean registration errors were obtained. The green dot shows the theoretical values of E and γ taken for the computation of results in Tables II and III.



  T e G Te D e G e R e , where G e is the strain-displacement matrix, D e is the stressstrain matrix and R e is a rotation matrix from corrotational formulation. The matrix D e is defined asD e =where λ and µ are the Lamé coefficients. The matrix G T e is defined as G e = ...

  Non-rigid Contour-based Registration of Cell Nuclei in 2D Live Cell Microscopy Images Using Dynamic Elasticity Model Dmitry V. Sorokin, Igor Peterlik, Marco Tektonidis, Karl Rohr, and Pavel Matula

TABLE III REGISTRATION

 III ERROR FOR LINE FEATURES, INNER POINTS AND END POINTS FOR DIFFERENT APPROACHES APPLIED FOR IMAGE SEQUENCES FROM THE SECOND EVALUATION DATASET.

						Mean error ēmean [pixels]
			Line features			Inner points			End points
	Sequence	B1	B2	B3	B4	B1	B2	B3	B4	B1	B2	B3	B4
	Unregistered	22.91 11.54 12.43 11.92	11.30 5.71 8.00 6.79 18.15 6.40 9.30 8.13
	Cont.-based [23]	24.35	9.13	9.02	11.58	14.60 5.50 5.33 4.57 15.94 3.44 4.92 7.63
	Int.-based [9]	25.99 11.16	9.28	9.99	22.47 7.50 6.56 4.43	9.42	4.56 5.01 2.51
	Cont.-based [25]	12.79	8.17	5.83	7.25	10.69 5.45 3.17 3.46	7.93	2.60 2.98 2.01
	Proposed	9.72	7.82	5.32	6.93	7.92	5.22 2.76 3.04	3.69	2.02 2.17

ACKNOWLEDGMENT

This work was supported by the Czech Science Foundation project 302/12/G157 and (ADD YOURS) .