Tobias Lieske
email: tobias.lieske@fau.de

Wilfried Uhring
email: wilfried.uhring@icube.unistra.fr

Norbert Dumas
email: norbert.dumas@icube.unistra.fr

Jérémie Léonard
email: jeremie.leonard@ipcms.unistra.fr

Anastasia Ioanna Skilitski

Dietmar Fey
email: dietmar.fey@fau.de

Embedded Fluorescence Lifetime Determination for High

Keywords: Throughput, Low-Photon-Number Applications FPGA, time-resolved fluorescence, fluorescence lifetime, microfluidics, embedded signal processing

come   L'archive ouverte pluridisciplinaire

Introduction

Fluorescence describes the emission of light from an electronically excited state. The lifetime τ is the average time between its excitation and the emission of a photon. Fluorescence spectroscopy and time-resolved fluorescence analysis methods are among the primary re-search tools in biochemistry and biophysics, but are not limited to these fields. The main biotechnological applications include flow cytometry, medical diagnostics, DNA sequencing, forensics and genetic analysis [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. A typical example that benefits from increased throughput is fluorescence-activated droplet sorting (FADS) [START_REF] Baret | Fluorescence-activated droplet sorting (fads): efficient microfluidic cell sorting based on enzymatic activity[END_REF] that is used for drug discovery [START_REF] Schmid | Sorting drops and cells with acoustics: acoustic microfluidic fluorescenceactivated cell sorter[END_REF].

Fluorescence spectroscopy measures the average emission intensity, but results underlie many influences such as turbidity, sample concentration, geometry of the sample and chromophore bleaching. Furthermore, this analysis method cannot capture important information, such as fluorescence emitted by one fluorophore masking other fluorophores due to spectral overlap as well as interactions of the fluorophore with the substrate or other macromolecules [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF].

Time-resolved fluorescence analysis methods on the other hand offer more information and are more robust, but require sophisticated optics and electronics due to the short time scale of fluorescence processes, which typically have emission rates of 10 8 s -1 [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. The analysis is based on the fluorescence lifetime, which has no dependency on emission intensity [START_REF] Giles | Characterization of erbium-doped fibers and application to modeling 980-nm and 1480-nm pumped amplifiers[END_REF][START_REF] Paschotta | Ytterbium-doped fiber amplifiers[END_REF] and is independent of the probe concentration as long as the probe emission is well above the background signal [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. The fluorescence lifetime enables a distinct classification.

Microfluidics is a technology that allows the processing of small amounts of fluids (10 -9 to 10 -18 liters), utilizing channels with diameters of tens to hundreds of micrometers. This analysis technique offers various benefits such as carrying out separations and detections with high resolution and sensitivity, low cost and short times for analysis [START_REF] Whitesides | The origins and the future of microfluidics[END_REF]. In order to combine the robustness of time-resolved fluorescence analysis with high throughput microfluidics, a sufficiently fast data processing is required to form a reliable high throughput sensing method. A particularly challenging example is microdroplet or cell sorting. The proposed experiment setup is illustrated in Fig. 1. A low cost high repetition rate picosecond laser diode generator is focused with lens L 1 on the microdroplets in the channel, causing fluorescence. The emitted photons (dotted lines) pass the dichroic mirror unimpeded, are collected by lens L 2 and then hit the single photon avalanche diode (SPAD). A time to digital converter (TDC) measures the time between the excitation and the emission of the photon. Based on this information, a histogram of arrival times can be generated, see Fig. 2. The fluorescence intensity is defined by Eq. 1, where I 0 is the initial intensity and τ the fluorescence lifetime [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. The histogram of arrival times in Fig. 2

I(t) = I 0 e -t τ (1)
lifetime is then used to separate the classified microdroplets correctly further back in the channel. Due to the flow velocity, the fluorescence lifetime computation has a hard real-time deadline set at 5 ms. The challenge of this approach is the complexity of the time-resolved fluorescence analysis method. Microfluidics systems have no problem generating a laminar high throughput flow of up to twenty thousand so called microdroplets per second [START_REF] Teh | Droplet microfluidics[END_REF]. However, there are currently two problems that limit the actual throughput of such a sorting system. The photon counting rate for current time-correlated single-photon counting (TCSPC) systems is limited by the sample's lifetime, the pileup effect and dead-times of electronics, which will all be discussed in Section 3.

In a previous experimental work, we have investigated how the pile-up effect and the photon shot noise -inherent to single photon detection -may limit the applications of TCSPC in very high throughput conditions [START_REF] Léonard | High-throughput time-correlated single photon counting[END_REF]. More specifically, we have demonstrated that these limitations do not preclude the use of TCSPC to detect with high sensitivity and good accuracy the fluorescence lifetime in microfluidic droplets of 50 microns diameter circulating at a flow rate exceeding 10 3 droplets per second. Recently, we applied the same experimental concept to demonstrate the feasibility and very good reliability of a biologically-relevant enzymatic activity assay based on fluorescence lifetime detection, in 200-picoliter droplets circulating at 300 droplets per second, by implementing TCSPC [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. In these previous experimental demonstrations, the data analysis (automatic droplet detection, building and fitting of fluorescence histograms for individual droplets) was done a posteriori, by post-processing the acquired data.

In this paper, we aim at enabling the real-time data processing in such experiments. This would be of central relevance to extend the above proof-of-principle experiments to microfluidic droplet sorting or fluorescence assisted cell sorting (FACS) experiments based on fluorescence lifetime detection by TCSPC. More specifically, we demonstrate the FPGA implementation of a monoexponential decay fit in order to achieve the processing of up to 151 histograms per millisecond, in conditions of low photon number per histogram.

Established methods for fluorescence lifetime computation require large photon counts for precise estimation, resulting in a low throughput. Algorithms based on a maximum likelihood estimator (MLE) deliver reliable high precision results even with low photon counts per histogram. This drastically reduces the required measurement time and in turn increases the possible input data rate. Transferring the data to a powerful workstation consumes too much time and embedded off-the-shelf products cannot handle the resulting data rate. Therefore, a dedicated low latency, high throughput embedded hardware architecture for highly accurate fluorescence lifetime classification based on low signal-to-noise inputs has been developed. This new architecture renders high throughput microfluidic sorting possible with current TCSPC systems. The system has been introduced in [START_REF] Lieske | Embedded fluorescence lifetime determination for high throughput real-time droplet sorting with microfluidics[END_REF] and is discussed here in more detail and validated with real experimental input data.

The following section reviews systems and methods for fluorescence lifetime extraction. After that, factors limiting the counting rate of TCSPC electronics for fluorescence measurements are explained in Section 3 and the MLE algorithm is revised in Section 4. Subsequently, the hardware architecture is explained in Section 5 and compared to existing approaches in Section 6. The paper ends with a validation of the presented hardware architecture with experimental input data in Section 7 and a brief conclusion in Section 8.

Related Work

The following review focuses on the fluorescence lifetime extraction from a given histogram of arrival times, as displayed in Fig. 2.

The work presented in [START_REF] Franch | A low cost fluorescence lifetime measurement system based on spad detectors and fpga processing[END_REF] shows a complete system similar to Fig. 1. While the TDC and histogram generation is implemented on a FPGA, the fluorescence lifetime is being computed externally on a workstation. Therefore, the histogram data is transferred via UART (Universal Asynchronous Receiver Transmitter). The results mention a delay of a few seconds to measure the fluorescence decay of a substance, which is clearly too slow for the real-time control system required by the sorting experiment proposed in this paper. Another approach using external processing is introduced in [START_REF] Gorpas | Real-time visualization of tissue surface biochemical features derived from fluorescence lifetime measurements[END_REF]. The fluorescence signals are transmitted via TCP/IP (Transmission Control Protocol/Internet Protocol) to a workstation that performs the lifetime computation. The signal acquisition, preprocessing and TCP/IP transmittance time is given as 9.5 ms. This time does not include the actual computation of the lifetime, proving that external computation of the fluorescence lifetime does not meet the given time constraints and an embedded approach is necessary to avoid long data transfer times.

There are several algorithms available that calculate the fluorescence lifetime. These algorithms can be classified into two groups. The first type are model based fitting methods such as Bayesian [START_REF] Rowley | Bayesian analysis of fluorescence lifetime imaging data[END_REF], maximum likelihood [START_REF] Maus | An experimental comparison of the maximum likelihood estimation and nonlinear least-squares fluorescence lifetime analysis of single molecules[END_REF] or maximum entropy [START_REF] Swaminathan | Analysis of fluorescence decay by the maximum entropy method: Influence of noise and analysis parameters on the width of the distribution of lifetimes[END_REF], method of moments [START_REF] Hall | Better estimates of exponential decay parameters[END_REF] and least squares deconvolution with Laguerre expansion [START_REF] Liu | A novel method for fast and robust estimation of fluorescence decay dynamics using constrained least-squares deconvolution with laguerre expansion[END_REF][START_REF] Gorpas | Real-time visualization of tissue surface biochemical features derived from fluorescence lifetime measurements[END_REF]. These fitting methods return precise results, but have a relatively high computational complexity and easily converge to a local minimum. Furthermore, each model is derived from application knowledge such as the number of lifetimes contained in the data and thus the fitting approach is limited to the selected data [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF].

The second algorithm class are nonfitting methods. Examples are the rotational invariance technique introduced in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF], the phasor analysis [START_REF] Digman | The phasor approach to fluorescence lifetime imaging analysis[END_REF], Prony's method [START_REF] Zhang | Prony's method for exponential lifetime estimations in fluorescence-based thermometers[END_REF], the integral equation method [START_REF] Li | Real-time fluorescence lifetime imaging system with a 32 × 32 0.13µm cmos low dark-count single-photon avalanche diode array[END_REF] and the center of mass method [START_REF] Li | Video-rate fluorescence lifetime imaging camera with cmos singlephoton avalanche diode arrays and high-speed imaging algorithm[END_REF][START_REF] Tyndall | A high-throughput time-resolved mini-silicon photomultiplier with embedded fluorescence lifetime estimation in 0.13 µm cmos[END_REF]. These algorithms offer a greatly reduced computational complexity when compared to the fitting based methods. Yet, some approaches are limited to monoexponential decays. An evaluation of these approaches is given in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF]. As will be shown in Section 6.2, the estimation quality of nonfitting methods is subpar and model based fitting methods offer superior results.

However, the main goal is to speed up the sorting process while maintaining the measurement precision. [START_REF] Maus | An experimental comparison of the maximum likelihood estimation and nonlinear least-squares fluorescence lifetime analysis of single molecules[END_REF] shows that the maximum likelihood estimation (MLE) performs better on low signal-to-noise fluorescence decays than the popular least squares method and a classification of monoexponential decays is possible with very few photons [START_REF] Köllner | How many photons are necessary for fluorescence-lifetime measurements?[END_REF]. In addition to that, the computational complexity for fluorescence lifetime estimation using the MLE algorithm is significantly lower than for the least squares algorithm [START_REF] Hall | Better estimates of exponential decay parameters[END_REF]. Therefore, the MLE algorithm has been selected for further analysis, since collecting less photons per droplet results in a shorter measurement time and in a higher experiment throughput.

Counting Rates of TCSPC Electronics

As already mentioned, the photon counting rate of TCSPC systems is limited by three factors. For one, TCSPC electronics exhibits a dead-time after a photon has been detected. Thus, photons arriving later at the SPAD can be masked by photons arriving early. This effect is called pulse pileup and deforms the recorded histogram in such a way that the measured fluorescence lifetime is reduced. This problem is usually addressed by following the 1% rule that states the excitation energy of one laser pulse should be so low that only one electron is emitted out of 100 excitation pulses. However, this 1% rule is being considered as overcautious. Counting rates at 10% decrease the lifetime by less than 1% and even a 30% count rate results only in a modest change [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Therefore, the technically maximum pulse repetition rate of the used laser can limit the photon counting rates. On the other hand, the maximum pulse repetition rate in order to determine the lifetime properly depends on the fluorescence decay times of the sample. The time between the pulses should be at least four times the longest lifetime in the sample [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Another technical limit is the actual duration of the dead-time of the TCSPC electronics. Effective counting rates further drop by 50% [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Considering an example with a lifetime of 5 ns, a maximum laser repetition rate of 50 MHz is allowed. Targeting a 30% count rate results in 15 Mcps (counts per second), which can be handled by currently available TCSPC products [2]. As will be shown in Section 6.2, a total photon count of 500 per histogram is enough for precise lifetime extraction with an average error below 0.5 % using the MLE algorithm. Provided one achieves a sufficiently fast fitting rate, up to 30000 microdroplets could be sorted per second. Our objective is therefore to implement a fitting routine that may process up to several 10 4 histograms per second.

Fluorescence Lifetime Estimation

The maximum likelihood estimator for monoexponential fluorescence decay signals has been derived in [START_REF] Hall | Better estimates of exponential decay parameters[END_REF] and is given in Eq. 2. With x = T /τ , T the time step of

1 + (e x -1) -1 -m(e mx -1) -1 = ν (2)
one histogram bin, m the number of histogram bins, ν as defined in Eq. 3 and N i the number of photons in his-

ν = m i iN i m i N i (3)
togram bin i. The term ν can be seen as the histogram bin index with the "average" number of photons. The left-hand side of Eq. 2 is strictly monotone decreasing and admits a unique solution as long as the source is decaying. Therefore, a solution can be easily obtained by applying the Newton-Raphson method, as displayed in Eqs. 4 to 6.

L(x) =1 + (e x -1) -1 -m(e mx -1) -1 -ν = 0 (4) x i+1 =x i - L(x i ) L (x i ) ( 5
) =x i - 1 + (e x -1) -1 -m(e mx -1) -1 -ν m 2 e mx (me mx -1) -2 -e x (e x -1) -2 (6)
Depending on the measured time window length, the total number of detected photons and the fluorescence lifetime, there is a visible noise floor, which can be seen in Fig. 2. As already mentioned, the photon count for the high throughput sorting application is relatively low. Therefore, higher index photon counts are relatively speaking stronger influenced by this noise floor, as these bins generally have a very low photon count. The noise floor deforms the exponential decay and directly influences ν and thus the estimated fluorescence lifetime value τ , see Eq. 3. Therefore, an iterative noise removal process called auto-windowing has been developed to acquire an improved approximation of the fluorescence lifetime. Depending on the estimated fluorescence lifetime, the input histogram is truncated to remove noisy bins and the fluorescence lifetime is recomputed, until a stable result is acquired.

This algorithm has been implemented in C on an ARM ® Cortex ® -A53 quad-core processor. Processing one histogram of 1024 bins takes 152.365 µs. Assuming no delays caused by concurrent memory accesses, a total of 26.25 histograms per millisecond can be processed with this embedded processor, which is not sufficient to handle the data rate of 30 histograms per millisecond of the numerical example given in Section 3. However, TCSPC electronics continuously improve enabling higher counting rates. Additionally, measuring a lifetime less than 5 ns allows higher laser repetition rates resulting in even higher counting rates, which in turn cannot be handled by the ARM ® Cortex ® -A53 quad-core processor.

Hardware Implementation

Therefore, a dedicated hardware implementation has been developed in VHDL to avoid compute-bound fluorescence measurement systems and maximize the sorting throughput. The presented fluorescence lifetime computation system is displayed in the block diagram in Fig. 3 connected to the histogram and control unit, as shown in the experiment setup in Fig. 1. The input and output interfaces use asynchronous FIFO (first-infirst-out) buffers for clock synchronization. The actual computation system consists of a preprocessing stage and a computation stage. The first part addresses the input data and performs some preprocessing to compute ν for Eq. 6. The second stage implements the Newton-Raphson iteration of the estimator function Eq. 6 to compute the fluorescence lifetime τ . This pipeline-like system structure enables concurrent processing of input data.

Preprocessing

The preprocessing stage is designed as a streaming architecture that accepts each input histogram bin by bin and computes the ν value, which is required to solve Eq. 6. As displayed in Eq. 3, ν depends on the number of histogram bins. However, the histogram size changes when the auto-windowing approach truncates the histogram. Therefore, different ν i values will be required to compute the fluorescence lifetime with this iterative approach. In order to avoid any delays by on-demand computations, ν i is computed on the fly for every histogram bin index i and stored in a buffering system.

The computations in this stage are performed on fixed-point data that scales throughout the module depending on the configurable input data width. This approach reduces the hardware resource requirements as opposed to a floating-point based implementation. Furthermore, the critical path of the preprocessing stage is reduced. In combination with the separate clock regions for the preprocessing and computation stages, this increases the maximum input data frequency and therefore the system throughput. The measurements performed in Sections 6 and 7 show that the fixed-point computation does not impair a precise fluorescence life-time computation when compared to a floating-point implementation.

The preprocessing module in Fig. 3 is required to extract the actual histogram that represents the exponential decay function from the noisy input provided by the histogram generation module, see Fig. 2. The two main noise components in such a measured histogram are the noise floor caused by background light and spikes in the histogram caused by Poisson noise, which is inherent for photon counting processes [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. In order to compensate both noise sources, the start of the histogram is computed as the mean index position of the maximum, which is determined within a configurable window at the beginning of the input histogram. In the course of an experiment, temperature changes influence the TDC behavior and can thus gradually shift the beginning of the histogram. Therefore, the mean maximum index is only computed over a dynamically configurable window size of the last histograms to account for this behavior.

While the preprocessing can be implemented in a fully pipelined streaming fashion, both the autowindowing approach and the Newton-Raphson method are iterative algorithms. In order to maximize the throughput, the computed ν i values are stored in a buffer system as single-precision floating point values. The ν i values of different input histograms are stored in separate buffers slots and accessible to the second stage to enable concurrent computations of Eq. 6.

Fluorescence Lifetime Computation

The computation of the fluorescence lifetime is performed in the computation stage of the system, see Fig. 3. Due to the large dynamic range of the exponential terms in Eq. 6, floating-point values are used instead of fixed point. In order to avoid that the critical path of floating-point arithmetic degrades the performance of the preprocessing stage, separate clocks are used for each stage and clock synchronization is implemented in the ν i buffer system module.

The second stage is fully pipelined, although autowindowing and Newton-Raphson are both iterative algorithms. The ν i buffer system in Fig. 3 signals the auto-windowing module when a new complete set of ν i values is available. The auto-windowing module in turn creates a new computation job and dispatches it to the Newton-Raphson module. A computation job consists of a number of job properties that are listed in Table 1. These are the current iteration indices j and k for the auto-windowing and Newton-Raphson iterations, the current histogram size m j , which decreases with each auto-windowing iteration due to the truncation process, and the current ν value which depends on the current histogram size. Further runtime parameters are set experiment-wide when configuring the computation system. These are listed in Table 2. The iteration limits j limit and k limit can be set within the range of statically configured values J limit and K limit set at synthesis time. Apart from the histogram data itself, Eq. 4 requires information about the time step T of one histogram bin during the measurement and an initial value x 0 = T /τ 0 . Therefore, the time step T and an initial fluorescence lifetime τ 0 are part of the configuration as well. Nevertheless, the computed fluorescence lifetime empirically converges within ten Newton-Raphson iterations, regardless of τ 0 .

With the job properties listed in Tables 1 and2, the Newton-Raphson module displayed in Fig. 3 computes k limit iterations of Eq. 6 in order to solve Eq. 4 and get a first estimate of the fluorescence lifetime τ . Once the Newton-Raphson iteration reaches k limit , the computed fluorescence lifetime τ is passed to the auto-windowing module where the histogram size m j for the next autowindowing iteration is derived and the corresponding ν mj is loaded from the ν i buffer system to update the job state in Table 1.

Figure 3 shows that the Newton-Raphson module controls the NR-iteration pipe block, which is a fully pipelined implementation of Eq. 6. In order to avoid numeric instability, critical parts of the computation are performed in double-precision. By pipelining the NR-iteration pipe block, up to 130 fluorescence lifetime computations can be performed in parallel each identified by a tag to ensure unambiguous assignment of input histogram to output lifetime. In reality, the number of parallel computation jobs is bound to the statically configurable number of buffer slots in the ν i buffer system module that is required to synchronize the input histogram stream with the computation per- formance. This means, the fluorescence lifetime computation for preprocessed histogram data stored as ν i values in buffer slot l has to be completed by the computation stage before the preprocessing stage overwrites the ν i values in slot l with data from a new input histogram. Using the input bin data frequency and the worst-case iteration counts J limit and K limit , the minimum number of buffer slots can be computed easily, as will be shown in Section 6.

Using this system structure, multiple preprocessing stages can be connected to share one fluorescence lifetime computation stage to save resources.

Evaluation

Implementation Results

In order to evaluate the presented fluorescence lifetime computation system, a parameter set for the configuration has to be selected. Table 3 shows an excerpt of the static configuration parameters and their values used for the implementation results in Table 4 obtained by Quartus Prime 15.1.0. The Data Width defines the input bit width of a histogram bin represented as fixed point number. m max configures buffer sizes and sets the maximum supported histogram size. AW Iter max configures the counter bit width for the auto-windowing iterations. The actual iteration count can be set dynamically between 1 and AW Iter max . The same applies for the Newton-Raphson iteration count that is configured via NR Iter max . The #Buffer parameter sets up the ν i buffer system.

An FPGA-based TDC and histogram generation module is introduced in [START_REF] Franch | A low cost fluorescence lifetime measurement system based on spad detectors and fpga processing[END_REF]. Such a module can be implemented along with the here presented fluorescence lifetime computation system on the same FPGA chip, given the unused resources in Table 4. Input histograms are streamed to the preprocessing stage bin by bin sequentially. Thus, the reported maximum clock frequency in Table 4 for the preprocessing stage corresponds to a potential input data stream of 151 histograms per millisecond, each histogram 1024 bins wide. Such a histogram rate may occur when the fluorescence lifetime needs to be measured at several points along the microfluidics channel or when counting rates of TCSPC electronics further improve.

Filling 1 of the 16 buffers in the ν i buffer system takes about 6.59 µs (1024 cycles). Assuming a continuous input data stream, the time window to process the data just stored in this buffer is about 98.98 µs (15360 cycles), which is the time to fill 15 buffers. Running the computation with the highest iteration count of 10 auto-windowing iterations, each consisting of 10 Newton-Raphson iterations, takes 93.86 µs (13840 cycles).

This means that even though an iterative algorithm has been chosen for the fluorescence lifetime computation, the presented architecture is IO bound due to smart parallel processing strategies. The given configuration can easily process 151 histograms per millisecond with a maximum computation delay of about 100,46 µs per histogram on a low cost FPGA development board [START_REF]DE0[END_REF]. The throughput and BRAM resource utilization of the system can be scaled down, by reducing the size of the ν i buffer system. In case higher input data rates need to be processed, several preprocessing stages can be connected with some glue logic to one shared computation stage. lifetimes τ ∈ {1 ns, 3 ns, 5 ns, 15 ns} have been generated. The time step T and number of histograms bins m have been selected according to Table 5. Poisson noise has been added to the histograms, since this kind of noise is inherent for photon counting processes and also is the prominent kind of noise for low photon counts [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Further, the behavior for low photon counts needs to be analyzed. Therefore, histograms with a different total photon count N ∈ {500, 600, . . . , 1500} have been generated. Due to the Poisson noise, the total number of photons deviates and is not exact for each histogram. For every {τ, N } set, 1000 histograms have been generated and the corresponding lifetime has been computed with the MLE algorithm implemented in C and by simulating the VHDL design. Figure 4 lists the coefficient of variation (CV), which is the quotient of the standard deviation to the mean. This ratio is a standardized measure of dispersion. The results show that the noise introduced by the processing is in an acceptable range compared to the input noise (CV of Poisson distribution is 1/ √ N ) and decreases with increasing N . Slight variations are caused by the added Poisson noise. Some histograms are noisier than others and some sets include more of such "noisier" histograms. Similarly, N varies from histogram to histogram.

Quality of Computed Fluorescence Lifetimes

The non fitting algorithms analyzed in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF] are compared by their normalized bias (∆τ /τ ) and the F -value, which is the normalized precision defined as F = √ N • σ τ /µ τ . An ideal computation has F = 1, while realistic computations feature F > 1 and F 1 [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF]. The histograms used in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF] have m = 1024 and a time bin width T = 12.5 ns / 1024 ≈ 12.2 ps, which is a much finer temporal resolution than the one used for the test cases in this paper, see Table 5. Unfortunately, N is not given, but the number of photons in the first bin is, which is either 1000 when τ is being varied or τ is fixed at 3 ns and the number of photons in the first bin varies from 100 to 5000. In order to put these numbers into context, the mean photon count in the first bin where the exponential decay starts is listed in Table 6 for all test cases used in this paper. It is obvious that the test cases used in this paper feature a drastically lower total photon count N , which implies a much shorter measurement time per droplet and thus improves the throughput. Nonetheless, the resulting bias and F-value measurements presented in Figs. 5 and6 respectively, show that the MLE algorithm returns superior results.

The bias for all measurements plotted in Fig. 5 is clearly below 0.5%, while the bias reported in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF] for τ = 3 ns for the integral equation method (IEM) is about -1.6% and about 6.3% for the center of mass method (CMM). Both clearly exceeding the bias of the MLE. Further, the listed phasor analysis features a bias of about 0.8% (τ = 3 ns), while the presented fluorescence lifetime estimation via rotational invariance techniques (FLERIT) algorithm is very close to 0% bias (τ = 3 ns). Extracting an exact range from the plot is not possible due to the coarse scale. The CMM method is discarded for the further comparison considering its excessive lifetime dependent bias.

Looking at the F-values reported in [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF], it is clear that FLERIT, IEM and phasor exhibit a lifetime dependent precision. The authors list optimized regions for the analyzed algorithms, where F < 4. These re- gions are from 0.4 ns to 14 ns for FLERIT, 0.4 ns to 8.8 ns for IEM and 0.4 ns to 4.23 ns for phasor. A higher resolution plot for τ = 3 ns shows that the F-value for the two best algorithms FLERIT and IEM varies between 1.5 and 2.0 [START_REF] Yu | Estimation of fluorescence lifetimes via rotational invariance techniques[END_REF]. In contrast to this, the F-value measurements of the MLE algorithm given in Fig. 6 show no dependency on the lifetime and no test case exceeds a F-value of 1.22.

Application to Real Experimental Data

In order to assess the performances of the above fitting procedure and VHDL implementation with real data, we use experimental data recorded in a proofof-principle experiment demonstrating the feasibility of an enzymatic activity assay based on fluorescence lifetime detection, at very high-throughput in microfluidic droplets [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. More specifically, in this assay, a fluorescently labeled substrate is used such that if the enzyme is active, the fluorophore is cleaved from the substrate and its fluorescence lifetime increases. Two data sets named 1 and 2, were acquired by implement- Fig. 7: Exemplary histogram data acquired from the experimental measurement setup in [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF] with additional postprocessing to remove the non-zero offset.

ing TCSPC in microfluidic droplets (as in Fig. 1, see also [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]) at a throughput of 300 droplets per second, in conditions where about 2500 single photons were collected per droplet. Two successive experiments where performed to generate data set 1 and 2 in which the droplet contents were the bare substrate (data set 1 ) or the substrate cleaved by the enzyme (data set 2 ). A histogram of the time-stamped single photon detection events was constructed for each droplet, by off-line data post-processing. Figures 7a and7b display exemplary three histograms from data set 1 and data set 2 respectively. Here we present the analysis of the two sets of histograms with the presented VHDL fitting scheme. (a) Fluorescence lifetime distribution for data sets 1 and 2. Fig. 8: Distribution of the fluorescence lifetimes that were computed by the presented VHDL implementation. The input data are the measurements from [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF] with additional postprocessing to remove the non-zero offset.

The original experimental data featured a non-zero offset that results in an overestimation of the fluorescence lifetime with the presented MLE algorithm and thus had to be removed before processing the histogram data. Further improvement will be to include an automated procedure for offset subtraction online during data acquisition.

For easier comparison with the results analysis from [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF], the computed fluorescence lifetime values, which are computed with the same system setup described in Section 6, are plotted as histograms in Fig. 8a. Similar to the evaluation in Section 6, the parameters µ τ , σ τ , the coefficient of variation and the F-value are listed in Table 7, along with the reference values published in [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. The results for σ τ , as well as the related coefficient of variation and the F-value, are in good agreement with the reported values from [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. The matching of these parameters indicates a high precision of the fluorescence lifetime estimate by the presented computation system. However, the computed fluorescence lifetimes do not accurately match. The main reason for this is that in [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF], the fluorescence lifetime of each droplet was computed offline with a 5-parameter fit implemented via a nonlinear least-square minimizing Matlab routine to fit each individual histogram, whereas here we perform a singleparameter fit based on the above MLE algorithm. In principle, both fitting procedure must give identical results in the case of purely monoexponential decays with zero offset, and strictly poissoninan noise distribution.

Here, the original experimental data do feature a nonzero offset, which strongly influences the fluorescence lifetime estimate if it is not taken into account. Hence, for the present validation, we had to estimate and subtract this average offset before processing the histogram data, which reduces the discrepancy between the average lifetimes to less than 20% between both analysis schemes, as shown in Tables 7 and8. In addition, the fluorescence decays are not rigorously monoexponential as can be readily seen for instance in Figure 7a (see residual systematic glitch at 15 ns, induced by a slight defect in the FPGA-based TDC). We conclude that the average fluorescence lifetime extracted by both fitting procedure is not affected in the same way by the deviations from purely monoexponential decays.

Future improvements to enhance the accuracy of online fluorescence lifetime determination by the present MLE routine will be to include an automated procedure for offset subtraction during data acquisition, and to correct the residuals defects of the TDC.

Most importantly, while we just discussed the origins of inaccuracy (i.e. bias in the average lifetimes) we now underline the very good precision (width of the lifetime distributions) of the present method, which is similar to that obtained by the off-line 5-parameter fit of ref [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. This very good precision is certainly the most critical and important result the present work, even in the case of non perfect data sets, because the precision is the most important criterion in screening applications where the goal is to be able to distinguish between two distinct values. Finally, following the example in [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF], the droplet histograms were summed up 10 by 10 to simulate a 10-fold increase of the effective droplets exposure time during the measurement and reduce the Poisson noise by a factor of √ 10. The computed fluorescence lifetimes for these two newly created data sets are plotted in Fig. 8b and the results of the statistical analysis are listed in Table 8 together with the reference values from [START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF]. As expected, the precision, estimated by the σ τ parameter, improves by a factor of ∼ √ 10. The measurement results in Tables 7 and8 show that, remarkably, the presented MLE-based fluorescence lifetime computation system still produces very high precision results, even if the accuracy may be affected by imperfect data quality.

Conclusion

Time-resolved fluorescence analysis is among the primary research tools in biochemistry and biophysics [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF]. Increasing the throughput of FADS systems improves the efficiency of drug discovery and can help finding new cures for diseases. This paper presents the implementation of a MLE based fluorescence measurement system. The results show that the algorithm offers superior fluorescence lifetime estimations, even for very low signal-to-noise inputs, i.e. low photon count. This enables shorter measurement times, which in turn increases the throughput of FADS systems, while also improving the measurement quality.

The MLE algorithm features indeed a higher computational complexity than non fitting algorithms. However, this paper shows that a flexible low cost FPGA implementation is capable of processing high data rates in real time.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Schematic setup of a microfluidics sorting system based on time-resolved fluorescence.

Fig. 3 :

 3 Fig. 3: System overview of the presented fluorescence lifetime computation unit.

Fig. 4 :

 4 Fig. 4: Coefficient of variation for the C and VHDL implementation, as well for the input histogram.

Fig. 5 :

 5 Fig. 5: Normalized bias for C and VHDL implementation. The optimum is a bias of 0.

FFig. 6 :

 6 Fig. 6: F-Value for C and VHDL implementation.

  Exemplary histograms from data set 2.

  Fluorescence lifetime distribution for data sets 1 and 2, where 10 individual droplet histograms were summed up to one input histogram.

  Tobias Lieske received his master degree in Information and Communication Technology at the Friedrich-Alexander-University (FAU) Erlangen-Nürnberg in 2015. Currently he is a research fellow at FAU at the chair of computer architecture. There he focuses on low-power system design, architecture generation and accelerator design. Wilfried Uhring received the M.Sc degree in microelectronic and the Master's Degree in engineering physic in 1999 and the Ph.D. degree in optoelectronics from University of Strasbourg, France, in 2002. Since 1999, he worked on the field of high speed imaging. He designed several ultrafast optical detection devices such as streak camera and gated intensified camera with sub-nanoseconds to picoseconds resolution. Since 2003, he extend his research activity in integrated ultrafast optoelectronic CMOS devices such as solid-state streak camera and System on a Chip. In 2013, he joined the Engineering science, computer science and imaging laboratory (ICube) laboratory, university of Strasbourg and CNRS in France, where he managed the Heterogeneous System and Microsystem team (SMH) of the ICube laboratory. Norbert Dumas received his Ph.D. on MEMS design and test from University Montpellier II in 2005. He then worked for 2 years at University of Lancaster (UK) as a postdoc in the framework of a Network of Excellence. He was a research engineer at LIRMM/University Montpellier II. He joined the University of Strasbourg as an associate professor in 2012 and a member of the ICube Laboratory "Heterogeneous Systems and Microsystems" team. His research is mainly focussed on integrated sensors, actuators and microfluidics for lab on a chip.After an experimental PhD Work (2000-03, Laboratoire Kastler Brossel, Paris) and a Post-Doctoral fellowship (2004, Universiteit van Amsterdam) in quantum physics on ultracold gases, Jérémie Léonard was appointed as a CNRS researcher in 2005 at IPCMS (Strasbourg, France). Since then he has been using UV-Vis femtosecond spectroscopy to study ultrafast photoreactions of organic molecules in condensed phase. Among other research interests, he contributes the application of time-resolved fluorescence to the investigation of biomolecular interactions with both fundamental and applied targets. Prof. Dr.-Ing. Dietmar Fey holds a diploma degree in Computer Science from Friedrich-Alexander-University (FAU) Erlangen-Nürnberg, Germany. In 1992 he received a Ph.D. from FAU with a work on an investigation about Using Optics in Computer Architectures. From 1994 to 1999 he researched at Friedrich-Schiller-University Jena where he made his habilitation. From 1999 to 2001 he worked as lecturer at University Siegen before he became a Professor for Computer Engineering at University Jena. Since 2009 he leads the Chair for Computer Architectureat Friedrich-Alexander-University Erlangen-Nürnberg (FAU). His research interests are in parallel computer architectures, parallel programming environments, parallel embedded systems, and memristive computing.

Table 1 :

 1 Computation job properties.

	Job Property Value
	AW-Iteration	j
	NR-Iteration	k
	hist size	m j
	ν	ν m j

Table 2 :

 2 Experiment-wide job properties.

	Experiment Property	Value
	AW-Limit	j limit
	NR-Limit	k limit
	Time Step	T
	Initial Lifetime	τ 0

Table 3 :

 3 Configuration parameters for implementation results in

Table 4 .

 4 Data Width m max J limit K limit #Buffer

	32	1024	10	10	16

Table 4 :

 4 Implementation results on Altera Cyclone V 5CSEMA4U23C6.

	f stage 1	f stage 2	FF	ALM	DSP	BRAM
	155.18 MHz 147.45 MHz 29891 (47 %) 11937 (75 %) 71 (85 %) 228 (84 %)

Table 5 :

 5 Time step and histogram size for each τ test case.

	τ [ns]	1	3	5	15
	T [ps]	45	90	90	90
	m	512 512 512 1024

Table 6 :

 6 Mean photon count in the first bin of the exponential decay. .00 31.61 34.88 40.09 43.34 48.86 52.15 57.67 61.78 67.17 3 13.34 16.56 19.19 21.85 24.50 27.10 29.24 33.07 35.16 37.88 41.23 5 8.83 10.35 11.89 13.80 15.63 16.89 18.84 20.77 22.18 24.05 25.63 15 2.97 3.46 4.10 4.60 5.24 5.89 6.57 7.31 7.68 8.24 8.86

	τ [ns]	500 600 700 800 900 1000 1100 1200 1300 1400 1500
	1 22.38 27

Table 7 :

 7 Statistical properties for the computation results for data sets 1 and 2 from Fig.8aand the corresponding reference results from[START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF].

				VHDL			reference
		u [ns] sig [ns] CV [%] F-Value	u [ns] sig [ns] CV [%]
	data set 1	8.3	0.25	3.0	1.5	6.9	0.25	∼3.6
	data set 2	12.2	0.47	3.8	1.8	11.4	0.42	∼3.6

Table 8 :

 8 Statistical properties for the computation results for data sets 1 and 2 from Fig.8band the corresponding reference results from[START_REF] Ioanna Skilitsi | Towards sensitive, high-throughput, biomolecular assays based on fluorescence lifetime[END_REF].

										VHDL	reference
								u [ns] sig [ns] CV [%] F-Value	u [ns] sig [ns] CV [%]
					data set 1		8.1	0.07	0.9	1.4	6.9	0.08	∼1.2
					data set 2		12.1	0.14	1.1	1.7	11.2	0.13	∼1.2
		300	data set 1
		250	data set 2
	# sample	100 150 200
		50
		0	7	8	9	10	11	12	13	14
					fluorescence lifetime [ns]

This work was financially supported by the Research Training Group 1773 "Heterogeneous Image Systems", funded by the German Research Foundation (DFG).


Fitsum Abadi

Olivier Gimenez

Hans Jakober

Wolfgang Stauber

Raphaël Arlettaz

Michael Schaub
email: michael.schaub@vogelwarte.ch

Estimating the strength of density dependence in the presence of observation errors using integrated population models

Keywords: Bayesian, demographic parameters, density-dependence, integrated population model, Lanius collurio, observation error, population growth rate, population size, simulation

1. Assessing the strength of density-dependence is crucial to understand population dynamics, but its estimation is difficult. Because estimates of population size and demographic parameters usually include errors due to imperfect detection, estimations of the strength of density-dependence will be biased if obtained with conventional methods which in addition lack statistical power to detect density-dependence.

2. We propose a Bayesian integrated population model that combines different sources of demographic information (capture-recapture, population counts, and data on reproductive success) into a single model to study density-dependence. The model allows assessing the effect of density both on the population growth rate and on the demographic parameters while accounting for imperfect detection. We studied the performance of this model by simulation and illustrate its use with a case study on red-backed shrikes Lanius collurio.

3. Our simulation results showed that the strength of density-dependence is identifiable and estimated with good precision. The strength of density-dependence was estimated with higher precision when the integrated population model was used than when a conventional regression model, which ignores the observation error, was applied. As expected, the conventional regression model tended to overestimate density-dependence at the population and to underestimate density-dependence at the demographic level. The analysis of the red-backed shrike data revealed negative density-dependence at the population level most likely mediated by a density-dependent decline in adult survival.

4. This work highlights the potential of integrated population models in assessing densitydependence and its practical application in population studies.

Introduction

Factors that regulate population size are broadly classified into density-independent (stochastic) and density-dependent (deterministic) factors (Williams, Nichols & Conroy 2002, page 136).

Density-independent factors (e.g. weather) affect all individuals in a population in the same way regardless of population density. In contrast, the effects of density-dependent factors intensify as the population density increases. Depending on their effect, they can be differentiated as negative and positive density-dependent processes. Negative effects of density-dependence occur if a demographic rate decreases as density increases, whereas positive effects occur when both increase [START_REF] Sinclair | Density dependence, stochasticity, compensation and predator regulation[END_REF]. Intraspecific competition for resources (e.g. food, nest site) and predation are the main drivers of negative density-dependence [START_REF] Newton | Population Limitation in Birds[END_REF][START_REF] Gunnarsson | Experimental evidence for density-dependent survival in mallard (Anas platyrhynchos) ducklings[END_REF], whereas a low chance to find mates at low population density is a typical mechanism resulting in positive, or inverse, density-dependence (Allee effect) [START_REF] Saether | Life history variation, population processes and priorities in species conservation: towards a reunion of research paradigms[END_REF][START_REF] Courchamp | Inverse density dependence and the Allee effect[END_REF][START_REF] Morris | Measuring the Allee effect: Positive density dependence in small mammals[END_REF]. The two mechanisms can operate within the same population at different periods of time, depending on current density [START_REF] Courchamp | Inverse density dependence and the Allee effect[END_REF].

In order to understand better how populations are regulated, the effect of density on different age or stage classes should be quantified. Knowledge of density-dependence is also crucial in practical applications such as conservation or harvest regulations [START_REF] Hanski | Random walks in a metapopulation: How much density dependence is necessary for long-term persistence[END_REF][START_REF] Sinclair | Density dependence, stochasticity, compensation and predator regulation[END_REF][START_REF] Drake | Density-dependent demographic variation determines extinction rate of experimental populations[END_REF]. Many studies have used long-term time series data to estimate the strength of density-dependence [START_REF] Dennis | Density-dependence in time-series observations of natural populations: estimation and testing[END_REF][START_REF] Lande | Estimating density dependence in time-series of agestructured populations[END_REF]. The main principle is to study whether there is a negative relationship between population growth rate and population size. An important limitation of this approach is that it needs the assumption of a population census, i.e., a complete enumeration without errors in the counts. If an observation error is present and not accounted for, it inflates the type I error and misleadingly indicates the presence of density-dependence [START_REF] Shenk | Sampling variance effects on detecting density dependence from temporal trends in natural populations[END_REF][START_REF] Freckleton | Census error and the detection of density dependence[END_REF][START_REF] Knape | Estimability of density dependence in models of time series data[END_REF][START_REF] Lebreton | Assessing density dependence: Where are we left? Modeling Demographic Processes in Marked Populations[END_REF]. More recently, state-space models have been used to reduce or remove bias in parameters or functional forms of density relationships resulting from observation error [START_REF] De Valpine | Fitting population models incorporating process noise and observation error[END_REF][START_REF] Jamieson | Density dependence in North American ducks[END_REF][START_REF] Dennis | Estimating density dependence, process noise, and observation error[END_REF]. However, estimates are unbiased only if the observation error is relatively small [START_REF] Knape | Estimability of density dependence in models of time series data[END_REF]. A further limitation is that the analysis of population counts only reveals the effect of density at the population level and, consequently, the demographic mechanisms causing densitydependence remain unknown. Some studies have used long-term data to assess the impact of density on demographic parameters by regressing estimates of demographic parameters on population counts [START_REF] Coulson | The relative roles of density and climatic variation on population dynamics and fecundity rates in three contrasting ungulate species[END_REF][START_REF] Saether | Population dynamical consequences of climate change for a small temperate songbird[END_REF][START_REF] Barker | Measuring density dependence in survival from mark-recapture data[END_REF][START_REF] Paradis | Exploring density-dependent relationships in demographic parameters in populations of birds at large spatial scale[END_REF][START_REF] Barbraud | Climate and density shape population dynamics of a marine top predator[END_REF][START_REF] Tavecchia | Density-dependent parameters and demographic equilibrium in open populations[END_REF]. To get unbiased estimates of the strength of density-dependence this approach requires estimates of demographic rates and population counts that are not subject to observation errors. If there are observation errors (either on the demographic rates and/or on the population counts), the test for densitydependence will suffer from a lack of statistical power [START_REF] Lebreton | Assessing density dependence: Where are we left? Modeling Demographic Processes in Marked Populations[END_REF]. Recently, Schofield, [START_REF] Schofield | Flexible hierarchical mark-recapture modeling for open populations using WinBUGS[END_REF] applied the Jolly-Seber model to mark-recapture data to estimate the strength of density-dependence on survival. Using this model, estimates of survival and of population sizes are obtained while accounting for imperfect detection, and thus the power to detect density-dependence is maximized.

Despite the development of different methods for studying density-dependence, there are currently no approaches to assess the effect of density on all demographic parameters simultaneously while accounting for observation error. This is a drawback because the full demographic mechanisms of density-dependent population regulation cannot rigorously be studied.

Here we propose the use of a Bayesian integrated population model to study densitydependence in a coherent fashion. Integrated population models combine different sources of demographic data into a single model and provide estimates of all demographic parameters and population size with improved precision and that are free of observation error [START_REF] Besbeas | Integrating mark-recapturerecovery and census data to estimate animal abundance and demographic parameters[END_REF][START_REF] Brooks | A Bayesian approach to combining animal abundance and demographic data[END_REF][START_REF] Schaub | Use of integrated modeling to enhance estimates of population dynamics obtained from limited data[END_REF]Abadi et al. 2010a). Parameters that are unidentifiable from a single source can be estimated with these models [e.g. fecundity [START_REF] Besbeas | Integrating mark-recapturerecovery and census data to estimate animal abundance and demographic parameters[END_REF]; immigration rate (Abadi et al. 2010b)]. These models also are very appealing and powerful when the sample size is small (Abadi et al. 2010a). So far, integrated population models have been successfully used to estimate key demographic parameters from single-site demographic data [START_REF] Brooks | A Bayesian approach to combining animal abundance and demographic data[END_REF][START_REF] Schaub | Use of integrated modeling to enhance estimates of population dynamics obtained from limited data[END_REF][START_REF] King | Identifying and diagnosing population declines: a Bayesian assessment of lapwings in the UK[END_REF]Véran & Lebreton 2008;[START_REF] Reynolds | Integrated data analysis in the presence of emigration and mark loss[END_REF]Abadi et al. 2010b) as well as from multi-site data [START_REF] Borysiewicz | An integrated analysis of multisite recruitment, mark-recapture-recovery and multisite census data[END_REF][START_REF] Péron | Studying dispersal at the landscape scale: efficient combination of population surveys and capture-recapture data[END_REF][START_REF] Mccrea | Multi-Site Integrated Population Modelling[END_REF] and their performance has been investigated using simulations [START_REF] Besbeas | Completing the ecological Jigsaw. Modeling Demographic Processes in Marked Populations[END_REF]Abadi et al. 2010a).

The model we propose for studying density-dependence potentially has several advantages. Firstly, because integrated population models involve the use of state-space models for the population counts, it is possible to assess density-dependence based on estimates of population indices corrected for random observation errors rather than on the population counts including errors, which avoids the confounding effect of observation errors. Secondly, because demographic parameters (e.g. age-specific survival, fecundity) are estimated, it allows a test of densitydependence for them, and thus the demographic mechanisms of density-dependent population regulation can be identified. It also allows assessing density-dependence at the population level using the derived population growth rate. Thirdly, combined analysis of demographic data results in improved precision of estimates, which is expected to increase the statistical power to detect density-dependence.

We conducted simulations to examine the performance of integrated population model in estimating the strength of density-dependence under different magnitudes of observation error. We illustrate the method using a data set on a population of red-backed shrikes Lanius collurio.

Previous studies have shown that the studied shrike population is regulated by density at the level of the population [START_REF] Pasinelli | Impact of density and environmental factors on population fluctuations in a migratory passerine[END_REF], but that neither survival nor fecundity were a function of local density (Schaub, Jakober & Stauber in press). The latter study was performed with conventional methods which are expected to have lower power than the integrated analysis. Our objective here was to evaluate whether density-dependence at the level of the demographic rates could be identified using the new proposed model.

Methods

We present first the integrated population model to estimate the strength of density-dependence, and then describe the simulation set up and the procedure to generate the data sets. Finally, we describe the case study and how the model is implemented in the Bayesian framework. We were motivated by a short-lived bird species from which we sampled capture-recapture and reproductive success data, as well as annual population counts. However, the estimation of density-dependence using the integrated population model is not restricted to this life history type and these data sets, other life histories and data sets could also be modelled by adaptations of the underlying population and observation models.

Modelling density-dependence in the integrated population model

Integrated population models first require the formulation of the likelihood for the available demographic data sets. Detailed discussion of this likelihood formulation is provided in Abadi et al. (2010a). In summary, we used the standard Cormack-Jolly-Seber model (CJS; [START_REF] Lebreton | Modelling survival and testing biological hypotheses using marked animals -A unified approach with case-studies[END_REF] for the capture-recapture data. We assumed the number of individual encounter histories summarized in the m-array (m) follows a multinomial distribution with parameters that are functions of age-specific apparent survival ( jv  , ad  ) and recapture probabilities (p). jv  denotes the probability that a newborn individual survives and becomes a yearling (1 year old) and ad  denotes the survival probability of yearlings and adults (older than 1 year). We denoted the likelihood for this model by cr

L ( | , , )

jv ad  mp . For the data on reproductive output, we assumed that the number of offspring produced per female (J) is Poisson distributed whose parameter is the product of fecundity (f) and the number of females recorded to produce young (R). Fecundity (i.e.

the number of newborns per adult females) of both yearlings and adults is assumed to be the same.

The likelihood for this model is denoted by L rp (J, R  f).

To describe the population counts we use a state-space model, which consists of a process and an observation model [START_REF] Besbeas | Integrating mark-recapturerecovery and census data to estimate animal abundance and demographic parameters[END_REF]. The process model describes the evolution of the underlying population sizes over time, thus it determines the link between the demographic rates and population size. Let N 1,t denote the number of 1 year old individuals in year t, N 2+,t denote the number of individuals older than 1 year in year t, and N t = N 1,t + N 2+,t , denote the total population size in year t. We assume that N 1,t+1 is generated by a Poisson process, as N 1,t+1 only takes a value between 0 (if no recruitment occurs) and ∞ (if recruitment is massive). The parameter of the Poisson process is the product of fecundity, juvenile survival and population size in year t. The sex ratio of newborn is assumed to be even and therefore we divide f by 2, and thus,

, 1, 1 ) ~Po( ( ) 2 t t jv t t f NN   eqn 1.
We further assume that N 2+,t+1 is generated by a Binomial process, as N 2+,t+1 only takes a value between 0 (if no individual survives) and the total population size a year before (if all survive). Thus, 2 , 1 , ~( , )

t t ad t N Bin N   eqn 2.
The observation model links the population counts (i.e. the number of detected and counted female breeders, denoted by y t ) to the latent population sizes (N t ). We then assume

y t ~ N(N t , 2 y  ) eqn 3
where the variance 2 y  quantifies the observation error. The likelihood for the state-space model is given by the product of the likelihood of the process and observation models, sy L ( | , , ) A simulation study by Abadi et al. (2010a) that combined and analyzed completely dependent data of capture-recapture, population counts, and reproductive success showed that the violation of the independence assumption only had a negligible effect on the precision of parameter estimates. We assume that this is true also in the present study.

jv ad  Nf  L ob (y  N,  2
All demographic parameters and population sizes appear in the likelihood of the integrated models as shown in eqn 4. This offers the possibility to model density-dependence of the demographic rates within a single model. We modelled the recapture probability with a random year parameter and assessed the effect of density on both survival probabilities and fecundity, , , * , 2 01 logit( )= ~(0, )

jv jv t jv t jv t t N N           eqn 5 , , * 2 , 23 logit( ) = ~(0, ) ad ad ad tt ad t t N N           eqn 6 * 2 5 4 log( ) = ~(0, ) t t t t f f f N fN       eqn 7
where * t N is the standardized population size (i.e., * () () We calculated this regression outside the integrated population model, because we did not want to induce any constraints on the density-dependence at the population level, since densitydependence at the population level is always the result of density-dependent variation of demographic rates.

tt t t N mean N N sd N   ), t N = 1,t N + 2, t N  is the
We also estimated the strength of density-dependence, ignoring the observation error, by regressing demographic rates estimated with single data analysis (i.e. survival probabilities from capture-recapture data estimated with the CJS model, and fecundity estimated from data on reproductive success with a Poisson regression model) on standardised population counts. The relationships between demographic parameters and population counts were similar with eqns 5-8, but here we used the population counts * t y rather than * t N and the population growth rate was calculated as  obs,t = y t+1 /y t and modelled as a function of * t y .

Simulation study

Using simulations we evaluated the identifiability of the estimates of the strength of densitydependence () and assessed their precision at varying levels of observation error under the integrated population model developed above. We also evaluated the effect of ignoring the observation error on the estimates of density-dependence obtained by regressing demographic rates on population counts.

A. Individual-based simulation of the data

The way we conducted simulations is analogous to that described in Abadi et al. (2010a). It consists of creating a population by simulating fates of individuals, sampling demographic data from this population, and analyzing these data with the developed integrated population model.

These steps are repeated several times and point estimates stored, serving the base to evaluate the performance of the model.

To create the population, we mimic a short-lived bird species and the development of the population sizes from one year to another is described by eqns 1 and 2. We specified the relationship between demographic parameters and the actual population size in year t as If the population size (N) increases by 10 individuals, the slope coefficients can be interpreted as the expected change in juvenile survival odds ratio to be exp(-0.5); in adult survival odds ratio to be exp(-0.2); and in fecundity to be exp(-0.1), indicating a negative effect of density on demographic parameters. The simulation parameters in eqns 9-11 were chosen in such a way that they result in reasonable demographic rates for a short lived bird species at the population size of about 50.

Next, we created the life history of individuals in the population for 20 years. To start the creation of the population, we set the initial population size at 40 individuals. We then specified the actual probabilities of survival and fecundity applying eqns 9-11 under consideration of the actual population size. For each individual alive in the population it was then simulated whether it survived for another year, how many juveniles it produced, and how many of the juveniles survived. The sum of the number of yearlings and surviving adults gives then the population size in the following year, and we specified the demographic rates for this following year. This was repeated for 20 years. Information about the state of an individual (dead, alive, in which age group it is given alive) and the number of offspring it has produced were stored. To avoid transition effects of the initial conditions we only used the last 10 years of the simulated data sets.

Once we created the population, the next step was to sample individuals for the different studies. We independently selected 1000 individuals at random from the population to be available for capture-recapture and reproductive success data sampling. To create the capture-recapture histories and reproductive success data based on individuals subject to each study, we set the initial capture probability at 0.90 and 0.50 for juveniles and adults (1 year and older), respectively, and the recapture probability at 0.50. The probability to detect and record reproductive success was set at 0.90. All these capture/recapture and detection probabilities were assumed to be constant across time. To create the population counts in each year, we used a binomial distribution with parameters actual population size in the given year (N t ) and detection probability (P s ). The detection probability was assumed to be constant over time. We considered two scenarios, one with detection probabilities of female breeders of 0.50, and one with 0.90. The lower detection probability corresponds to a large observation error in the population counts. This is because the variance of binomial sampling (i.e. NP s (1-P s )) is largest with P s =0.50. Our sampling procedures did not result in completely independent data sets as some individuals that were involved in either capture-recapture or reproductive success data sampling also had a chance to be included in the population counts. We simulated 500 data sets and analyzed them using the proposed integrated population models.

B. Case study

We was about 50 pairs. For the analysis we assumed an even sex ratio and only considered capturerecapture data of females. The goal was to estimate the strength of density-dependence at the population level and of the demographic rates to assess 1) whether the population is regulated by density, and 2) if so, through which demographic mechanisms.

Because the studied red-backed shrike population is geographically open, we had to extend the integrated population model to include immigration (Abadi et al. 2010b). Therefore, we modified the state-space model (eqns 1-3) as follows: the number of 1 year old individuals N 1,t+1 followed a Poisson process with 1, 1 , ~Po( ( ) ) 2

t t t jv t f NN   eqn 12
The number of immigrants N im,t+1 was modelled with a Poisson distribution as ,1 ()

im t t t N Po N   : eqn 13
where the immigration rate t  is defined as the number of female immigrants of age 1 year or older in year t+1 per breeding females in the previous year N t .

The number of at least 2 years old surviving individuals followed a binomial process with , 1 , ~( , )

ad t t ad t N Bin N   eqn 14
The total number of individuals in year t is given by 1, , ,

t t ad t im t N N N N    .
For the observation equation, we assumed a normal distribution, hence y t the annual number of counted breeding females in year t distributed as a N(N t , 2 y  ) .

To estimate the strength of density-dependence on the demographic parameters as well as the population growth rate, we fitted the models specified in eqns 5-8. We also modelled immigration rate with a random year parameter, but did not impose density-dependence on it.

Further, the recapture probability (p) was modelled time-dependent where time was treated as random. Beside the model that includes density dependence for the demographic parameters, we also fitted a model without density-dependence (i.e. 1 3 5 0       in eqns 5-7), in order to get unconstrained annual estimates of the demographic parameters.

C. Bayesian analysis of the model Simulation study

Our simulation study clearly showed that the strength of density-dependence in the demographic rates as well as in the population growth rate was identifiable with the integrated population model regardless of the magnitudes of observation error we considered, indicated by smooth posterior densities with a clear peak (Figs 1 &2). If the parameters were not identifiable, we would have expected a posterior density that is similar to the specified prior distribution [START_REF] Gimenez | Weak identifiability in models for markrecapture-recovery data[END_REF]). The widths of the posterior distributions resulting from the integrated population model were narrower when the observation error was small compared to when it was large. This indicates that the precision of the density-dependence estimator declines with increasing observation error of the population counts. In general, the strength of density-dependence was estimated with more precision at the population level than at the demographic level, indicated by a narrow width of the posterior distributions (Figs 1 & 2).

Estimates of the strength of density-dependence obtained from the integrated population model were generally more precise compared to that of the regression analysis. This was particularly evident for adult survival (Figs 1 & 2). As expected, the estimates of strength of density-dependence in the demographic parameters were slightly shifted towards zero when the observation error was large and not accounted for. However, this effect was not very strong. The model which ignores the observation error yielded stronger and slightly less precise estimates of density-dependence at the population level in comparison to the integrated population model (Figs

1 & 2).

Case study

Since there is no established goodness-of-fit test of the complete integrated population model available, we tested the goodness-of-fit just for the capture-recapture model using contingency tables (Pradel et al. 2005) via program U-CARE [START_REF] Choquet | U-CARE: Utilities for performing goodness of fit tests and manipulating CApture-REcapture data[END_REF]): there was no sign of lack of fit ( 2 59  = 42.75, P = 0.94). Based on the integrated population model without density- dependence, the annual estimates of age-specific survival, fecundity and population growth rate tended to decline with increasing population size in red-backed shrikes (Fig. 3), suggesting density-dependence in these demographic parameters and at the population level. We then explicitly estimated the strength of density-dependence in the demographic parameters and in the population growth rate using the integrated population model with density-dependence. The estimated linear relationships between demographic rates as well as population growth rate (on the transformed scale) and population size are shown in Fig. 3. The posterior distributions of the strength of density-dependence were wide for juvenile survival and fecundity, and more pointed for adult survival and population growth. As expected, the posterior means of the strengths of density-dependence were negative in all parameters (Fig. 4). The posterior probability that the estimated effect of density were negative was high for population growth (0.960), followed by adult survival (0.792), fecundity (0.598) and juvenile survival (0.559). Thus, there was clear evidence of a density-dependent population regulation at the level of the population that was more likely due to adult survival, than due to juvenile survival or fecundity.

Discussion

Density-dependence is an important ecological concept and understanding how it operates is crucial in conservation, harvesting and for accurate demographic projections [START_REF] Sinclair | Density dependence, stochasticity, compensation and predator regulation[END_REF][START_REF] Lande | Estimating density dependence in time-series of agestructured populations[END_REF]. In this paper we provide a framework for studying density-dependence using a Bayesian integrated population model. A key advantage of this unifying framework is that it allows estimating strength of density-dependence both at the demographic and population level while accounting for observation error and, consequently, the demographic mechanism causing density-dependent population regulation can be identified. In contrast, most existing techniques for demographic rates is weak, all demographic rates work jointly, and thus the effect at the population level magnifies. Second, the red-backed shrike population is geographically open, and immigration is substantial (mean (sd): 0.545 (0.036)). It is well possible that immigration is regulated by density, i.e., few individuals immigrate in years where the number of survivors and local recruits is high, and vice-versa. Yet, immigration is a parameter in our integrated population model that is estimated without having explicit data (Abadi et al. 2010), and for such parameters it appears that the strength of density-dependence cannot be estimated (i.e., it is not an identifiable parameter of the model). To evaluate for which demographic parameter density-dependence is important, we computed the probability that the strength of density-dependence is negative and use this as a testing criterion. Alternatively, one could also apply model selection for the same purpose. The set of models could then include models that impose density-dependence in some demographic parameters only. Yet, model selection in Bayesian hierarchical models is not an easy task [START_REF] Link | Model weights and the foundations of multimodel inference[END_REF][START_REF] Millar | Comparison of hierarchical Bayesian models for overdispersed count data using DIC and Bayes' factors[END_REF]).

In our simulation as well as in the case study we expressed density in terms of population size. Since the study area associated with the red-backed shrike population remained the same over time, the use of population size as a measure of density was justified. However, it is important to note that the influence of density on biological processes can be overlooked due to an inappropriate measure of density (Williams et al. 2002;[START_REF] Barker | Measuring density dependence in survival from mark-recapture data[END_REF]. Moreover, the main driving force for density-dependence is often competition for resources and not space. Ideally, one would therefore model the strength of density-dependence not with population size, but with the available resources per individual. Our model could be extended in this way if an estimate of resource availability could be obtained.

The goals of our simulation were mainly on the assessment of the identifiability of the estimates of strength of density-dependence and on the evaluation of the effect of observation error. We specified a short study period (i.e., 10 years), thus conditions where the estimation of density-dependence were difficult. Yet, the model was able to estimate density-dependence with good precision in this set-up. We could further examine the performance of the model in different directions. For instance, we could set up a simulation varying the strength of density-dependence and also the length of the study period. The ability to detect density-dependence often increases with long time series data [START_REF] Brook | Strength of evidence for density dependence in abundance time series of 1198 species[END_REF]. In our simulation study, we focused on the common kind of observation error (i.e. non-detection or false-negative error). One might further evaluate the performance of the model considering other kinds of observation error such as falsepositive errors occur.

The integrated population model is very flexible to include different shapes of densitydependence. Here we used the density-dependence model of the Ricker type for population growth rate [START_REF] Dennis | Density-dependence in time-series observations of natural populations: estimation and testing[END_REF], but specifying other density-dependence models which are discussed in the literature (e.g. [START_REF] Dennis | Density-dependence in time-series observations of natural populations: estimation and testing[END_REF][START_REF] Jamieson | Density dependence in North American ducks[END_REF]) is straightforward.

Moreover, some studies have shown that the effect can be non-linear in the parameters [START_REF] Paradis | Exploring density-dependent relationships in demographic parameters in populations of birds at large spatial scale[END_REF][START_REF] Saether | Pattern of variation in avian population growth rates[END_REF][START_REF] Tavecchia | Density-dependent parameters and demographic equilibrium in open populations[END_REF]. The integrated population model could also be extended to assess non-linear density-dependence, for instance by using penalized splines [START_REF] Gimenez | Semiparametric regression in capture-recapture modelling[END_REF]). In the present study, we looked at the direct effect of density but it is possible that density-dependence can affect demographic parameters and population growth rate with time lag [START_REF] Paradis | Exploring density-dependent relationships in demographic parameters in populations of birds at large spatial scale[END_REF]. Our model can easily be extended to study delayed densitydependence.

Overall, our model provides an improved statistical tool over current methods for estimating density-dependence. We believe it has great potential in conservation, management and ecology.

  y ). The likelihood of the complete integrated population model is obtained as the product of the likelihoods of all three data sources under the assumption of independence,

  population size in year t , and the  2 are the temporal variances of the demographic parameters as well as of recapture probability. The goal is to estimate the regression coefficients (, ) and the magnitude of temporal variability ( 2 ). The slope coefficients ( 1 ,  3 ,  5 ) estimate the strength of density-dependence on juvenile survival, adult survival and fecundity, respectively. For the recapture probability, which is a nuisance parameter, we assumed a random year effect, the strength of density-dependence on the population growth rate (i.e. at the level of the population). This was done by first calculating the population growth rate ( t = N t+1 /N t ), and then regressing the growth rate on population size as

  used 26 years of demographic data of red-backed shrikes collected from 1981-2006 in southwestern Germany near Göppingen (4839'N, 947'E) to illustrate the method. The red-backed shrike is a medium-sized (~30g) migratory passerine living in agricultural landscapes during the breeding season (Glutz von Blotzheim & Bauer 1993). The study area has a size of about 18 km 2 and mostly consists of meadows and pastures with interspersed orchards and hedgerows (for a more detailed description see Jakober & Stauber 1987). During numerous visits we localised each year the nests of nearly all breeding pairs, recorded the number of fledglings of each nest (n = 3580 fledglings from 1242 nests) and ringed the nestlings at an age of about 8 days (n = 3598). Territorial adults were caught with clap or mist nets and ringed with individual colour rings (n = 513 females). Each spring we visited the study area almost daily to resight colour-marked individuals. The population size varied between 35 and 71 observed breeding pairs, the mean size

Fig. 4 .

 4 Fig. 4. Posterior distributions of the estimated strength of density-dependence on juvenile survival

Acknowledgments

We are grateful to M. Kéry, D. Fletcher and two anonymous reviewers for useful comments on the manuscript. We acknowledge the financial support of the Swiss National Science Foundation (grant No. A0-107539 to Michael Schaub) and SATW (Germaine de Staël).

All the models were fitted within the Bayesian framework, specifying non-informative priors in order to reflect limited knowledge about the parameters. Specifically, we assigned a N(0, 100) distribution on the regression coefficients ( and  ). A N(100, 100) distribution truncated to positive values was assigned for age specific initial population sizes, and a U(0, 10) distribution for the temporal standard deviations of demographic parameters. To assess convergence of the Markov chain Monte Carlo (MCMC) algorithm to the stationary distribution, we randomly chose a single simulated data set and run three chains, each with different initial values, of 20000 iterations with a burn-in of 15000 iterations. The R [START_REF] Brooks | General methods for monitoring convergence of iterative simulation[END_REF] values were less than 1.1 for all parameters, suggesting convergence. We therefore run a single chain of 30000 MCMC iterations with a burn-in of 20000 thinning every 10 th observation in order to compute the posterior summary statistics. Using this setting, the analysis of one data set took approximately 50 minutes (on a 3.3 GHz processor, 2 GB RAM PC). For the case study, we first run three chains of 20000 iterations with a burn-in of 10000 to check whether convergence was reached. Since convergence was obtained (all the R values < 1.02), we run a single chain of length 100000, discarded the first 50000 as burn-in and thinned every 10 th observation. The run time to analyse the red-backed shrike data set was approximately 5 days (on a 3.3 GHz processor, 2 GB RAM PC). The posterior summary statistics were then computed based on 5000 samples. We used the R software version 2.9.1 (R Development Core Team 2008) to simulate the data and the analyses were done using the WinBUGS software calling it from R through the package R2WinBUGS [START_REF] Sturtz | R2WinBUGS: A package for running WinBUGS from R[END_REF]. The R and WinBUGS codes for the density-dependence model of the red-backed shrike data are provided in Appendix S1 in Supporting Information.

Results

testing density-dependence rely on time series data and do not account for observation error, or focus either on the effect of density on population growth rate or on a single demographic parameter only.

Using simulations, we showed that integrated population models provide estimates of the effect of population size on demographic parameters as well as on population growth rate with good precision. The parameters estimating the strength of density-dependence were identifiable with our model regardless of the magnitudes of observation error we considered. As expected, the estimates were slightly more precise when the observation error was small than large. Our simulation study also highlight that ignoring observation error had differential impact on the estimators of the strength of density-dependence at the population than at the demographic level.

The conventional regression model tended to estimate stronger density-dependence at the population than at the demographic level. These results are in agreement with previous findings (e.g. [START_REF] Shenk | Sampling variance effects on detecting density dependence from temporal trends in natural populations[END_REF][START_REF] Freckleton | Census error and the detection of density dependence[END_REF][START_REF] Lebreton | Assessing density dependence: Where are we left? Modeling Demographic Processes in Marked Populations[END_REF].

We also illustrated the application of the method using demographic data from a geographically open red-backed shrike population. We found strong support of densitydependence for the population growth rate. Of the considered demographic parameters, only adult survival was likely to be regulated by density-dependence, while juvenile survival and fecundity were hardly impacted by breeding density. Using the same data analysed with regression analyses with population counts uncorrected for observation error, Schaub et al. (in press) did not find support of density-dependence for survival and fecundity. With the application of the integrated population model the power to detect density-dependence increases (see simulations), which is the reason why we now found support of density-dependence operating on adult survival. One may wonder how there can be strong density-dependence for population growth and only relatively weak density-dependence for the demographic rates. First, even if density-dependence for the Williams, B.K., Nichols, J.D. & Conroy, M.J. ( 2002) Analysis and management of animal populations. Academic Press, San Diego, USA.

Supporting Information

The following supporting information is available for this article online.