Pouzat_ENP2018_handout.org

Fluorescence Imaging Analysis: The Case of Calcium Transients.

Codes

Define koma-article

Get the figure from van Vliet et al chapter

We start by downloading the PDF file of the chapter from the author’s website:

echo http://homepage.tudelft.nl/e3q6n/publications/1998/WaS98LVFBea/WaS98LVFBea.pdf

#+RESULTS[0e276d911c7ff39ffe1124ba6359f3aa835b52f5]: url-vanVlietEtAl-1998

http://homepage.tudelft.nl/e3q6n/publications/1998/WaS98LVFBea/WaS98LVFBea.pdf

wget $url 

#+RESULTS[9c9c2f685ca31d3d3a48fa299f10de4a73865190]: download-vanVlietEtAl-1998

convert WaS98LVFBea.pdf[4] -crop 500x285+50+150\! +repage figs/vanVlietEtAl_1998_Fig4.png

#+RESULTS[68d18d159b75b91e28364a7e4cc1bb822ed44ac5]: vanVlietEtAl-1998-extract-figure-4

Introduction

The variability inherent to fluorescence imaging data

ADU counts (raw data) from Fura-2 excited at 340 nm. Each square corresponds to a pixel. 25.05 s of data are shown. Same scale on each sub-plot. Data recorded by Andreas Pippow (Kloppenburg Lab. Cologne University).

One of the central pixels of the previous figure.

What do we want?

Given the data set illustrated on the last two slides we might want to estimate parameters like:

• the peak amplitude
• the decay time constant(s)
• the baseline level
• the whole time course (strictly speaking, a function).

If we have a model linking the calcium dynamics—the time course of the free calcium concentration in the cell—to the fluorescence intensity like: \[\frac{\mathrm{d}Ca_t}{\mathrm{dt}} \left(1 + κ_F(Ca_t) + κ_E(Ca_t) \right) + \frac{j(Ca_t)}{v} = 0 \, , \] where $Ca_t$ stands for $[Ca²⁺]_free$ at time t, $v$ is the volume of the neurite—within which diffusion effects can be neglected—and \[j(Ca_t) ≡ γ (Ca_t - Ca_steady) \, ,\] is the model of calcium extrusion—$Ca_steady$ is the steady state $[Ca²⁺]_free$ — \[κ_F(Ca_t) ≡ \frac{F_total \, K_F}{(K_F + Ca_t)^2} \quad \mathrm{and} \quad κ_E(Ca_t) ≡ \frac{E_total \, K_E}{(K_E + Ca_t)^2} \, ,\] where $F$ stands for the fluorophore en $E$ for the endogenous buffer.

Assuming that the fluorophore (Fura) parameters: $F_total$ and $K_F$ have been calibrated, we might want to estimate:

• the extrusion parameter: $γ$
• the endogenous buffer parameters: $E_total$ and $K_E$

using an equation relating measured fluorescence to calcium: \[Ca_t = K_F \, \frac{S_t - S_min}{S_max - S_t} \, ,\] where $S_t$ is the fluorescence (signal) measured at time $t$, $S_min$ and $S_max$ are calibrated parameters corresponding respectively to the fluorescence in the absence of calcium and with saturating $[Ca²⁺]$ (for the fluorophore).

• The variability of our signal—meaning that under replication of our measurements under the exact same conditions we wont get the exact same signal—implies that our estimated parameters will also fluctuate upon replication.
• Formally our estimated parameters are modeled as random variables and it is not enough to summarize a random variable by a single number.
• If we cannot get the full distribution function for our estimated parameters, we want to give at least ranges within which the true value of the parameter should be found with a given probability.
• In other words: an analysis without confidence intervals is not an analysis, it is strictly speaking useless since it can’t be reproduced—if I say that my time constant is 25.76 ms the probability that upon replication I get again 25.76 is essentially 0; if I say that the actual time constant has a 0.95 probability to be in the interval [24,26.5], I can make a comparison with replications.

A proper handling of the “variability” matters

Let us consider a simple data generation model: \[Y_i ∼ \mathcal{P}(f_i)\, , \quad i=0,1,\ldots,K \; ,\] where $\mathcal{P}(f_i)$ stands for the Poisson distribution with parameter $f_i$ : \[\mathrm{Pr}\{Y_i = n\} = \frac{(f_i)^n}{n!} exp (-f_i)\, , \quad \mathrm{for} \quad n=0,1,2,\ldots \] and \[f_i = f(δ i| f_∞, Δ, β) = f_∞ + Δ \, exp (- β \, δ i)\; ,\] δ is a time step and $f_∞$, Δ and β are model parameters.

Data simulated according to the previous model. We are going to assume that $f_∞$ and $Δ$ are known and that $(t_1,y_1)$ and $(t_2,y_2)$ are given. We want to estimate $β$.

Two estimators

We are going to consider two estimators for $β$:

• The “classical” least square estimator: \[ ˜{β} = arg min ˜{L}(β) \; ,\] where \[ ˜{L}(β) = ∑_j \big( y_j - f(t_j \mid β) \big)^2 \; .\]
• The least square estimator applied to the square root of the data: \[\hat{β} = arg min \hat{L}(β) \; ,\] where \[ \hat{L}(β) = ∑_j \big( \sqrt{y_j} - \sqrt{f(t_j \mid β)} \big)^2 \; .\]

We perform an empirical study as follows:

• We simulate 100,000 experiments such that: \[ (Y_1,Y_2) ∼ \big(\mathcal{P}(f(0.3|β_0), \mathcal{P}(f(3|β_0)\big) \; ,\] with $β_0=1$.
• For each simulated pair, $(y_1,y_2)^[k]$ ($k=1,\ldots,10^5$), we minimize $˜{L}(β)$ and $\hat{L}(β)$ to obtain: $(˜{β}^[k],\hat{β}^[k])$.
• We build histograms for $˜{β}^[k]$ and $\hat{β}^[k]$ as density estimators of our estimators.

Both histograms are built with 100 bins. $\hat{β}$ is clearly better than $˜{β}$ since its variance is smaller. The derivation of the asymptotic (large sample) densities is given in Joucla et al (2010).

CCD camera noise

CCD basics

Source: L. van Vliet et col. (1998) Digital Fluorescence Imaging Using Cooled CCD Array Cameras (figure 3).

“Noise” sources in CCD

• The “Photon noise” or “shot noise” arises from the fact the measuring a fluorescence intensity, λ, implies counting photons—unless one changes the laws of Physics there is nothing one can do to eliminate this source of variability (improperly called “noise”)—: \[\mathrm{Pr}\{N=n\} = \frac{λ^n}{n!} exp -λ\, , \quad n \, = \, 0,1,\ldots\, , \quad λ > 0\; .\]
• The “thermal noise” arises from thermal agitation which “dumps” electrons in potential wells; this “noise” also follows a Poisson distribution but it can be made negligible by cooling down the camera.
• The “read out noise” arises from the conversion of the number of photo-electrons into an equivalent tension; it follows a normal distribution whose variance is independent of the mean (as long as reading is not done at too high a frequency).
• The “digitization noise” arises from the mapping of a continuous value, the tension, onto a grid; it is negligible as soon as more than 8 bit are used.

A simple CCD model

• We can easily obtain a simple CCD model taking into account the two main “noise” sources (photon and read-out).
• To get this model we are going the fact (a theorem) that when a large number of photon are detected, the Poisson distribution is well approximated by (converges in distribution to) a normal distribution with identical mean and variance: \[\mathrm{Pr}\{N=n\} = \frac{λ^n}{n!} exp -λ ≈ \mathcal{N}(λ,λ) \; .\]
• In other words: \[ N ≈ λ + \sqrt{λ} \, ε \; ,\] where $ε ∼ \mathcal{N}(0,1)$ (follows a standard normal distribution).
• A read-out noise is added next following a normal distribution with 0 mean and variance $σ_R^2$.
• We are therefore adding to the random variable $N$ a new independent random variable $R ∼ \mathcal{N}(0,σ_R^2)$ giving: \[M ≡ N+R ≈ λ + \sqrt{λ+σ_R^2} \, ε \; ,\] where the fact that the sum of two independent normal random variables is a normal random variable whose mean is the sum of the mean and whose variance is the sum of the variances has been used.
• Since the capacity of the photo-electron weels is finite (35000 for the camera used in the first slides) and since the number of photon-electrons will be digitized on 12 bit (4096 levels), a “gain” $G$ smaller than one must be applied if we want to represent faithfully (without saturation) an almost full well.
• We therefore get: \[Y ≡ G ⋅ M ≈ G \, λ + \sqrt{G^2 \, (λ+σ_R^2)} \, ε \; .\]

For completeness: Convergence in distribution of a Poisson toward a normal rv

For completeness: Convergence in distribution of a Poisson toward a normal rv

We use the moment-generating function and the following theorem (e.g. John Rice, 2007, Mathematical Statistics and Data Analysis, Chap. 5, Theorem A):

• If the moment-generating function of each element of the rv sequence $X_n$ is $m_n(t)$,
• if the moment-generating function of the rv $X$ is $m(t)$,
• if $m_n(t) → m(t)$ when $n → ∞$ for all $|t| ≤ b$ where $b > 0$
• then $X_n \xrightarrow{D} X$.

Lets show that: \[Y_n = \frac{X_n - n}{\sqrt{n}} \; , \] where $X_n$ follows a Poisson distribution with parameter $n$, converges in distribution towards $Z$ standard normal rv.

We have: \[m_n(t) ≡ \mathrm{E}\left[exp(Y_n t)\right] \; ,\] therefore: \[m_n(t) = ∑_k=0^∞ exp\left(\frac{k-n}{\sqrt{n}}t\right) \frac{n^k}{k!} exp(-n) \; ,\] \[m_n(t) = exp(-n) exp(-\sqrt{n}t) ∑_k=0^∞ \frac{\left(n exp\left(t/\sqrt{n}\right)\right)^k}{k!}\] \[m_n(t) = exp\left(-n - \sqrt{n} t+ n exp(t/\sqrt{n})\right)\] \[m_n(t) = exp\left(-n - \sqrt{n} t+ n ∑_k=0^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!}\right)\] \[m_n(t) = exp\left(-n - \sqrt{n} t+ n + \sqrt{n} t + \frac{t^2}{2} + n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!}\right)\] \[m_n(t) = exp\left( \frac{t^2}{2} + n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!}\right)\]

We must show: \[n ∑_k=3^∞\left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} →_{n → ∞} 0 \quad ∀\ |t| ≤ b, \quad \text{where} \quad b > 0\, ,\] since $exp(-t^2/2)$ is the moment-generating function of a standard normal rv. But \[\left| n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} \right| →_{n → ∞} 0 \quad ∀\ |t| ≤ b, \quad \text{where} \quad b > 0\,\] implies that since \[- \left|n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} \right| ≤ n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} ≤ \left| n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} \right| \, .\]

But for all $|t| ≤ b$ where $b > 0$ \begin{displaymath} \begin{array}{lcl} 0 ≤ \left| n ∑_k=3^∞ \left(\frac{t}{\sqrt{n}}\right)^k \frac{1}{k!} \right| & ≤ & n ∑_k=3^∞ \left(\frac{|t|}{\sqrt{n}}\right)^k \frac{1}{k!}
& ≤ & \frac{|t|^3}{\sqrt{n}} ∑_k=0^∞ \left(\frac{|t|}{\sqrt{n}}\right)^k \frac{1}{(k+3)!} \ & ≤ & \frac{|t|^3}{\sqrt{n}} ∑_k=0^∞ \left(\frac{|t|}{\sqrt{n}}\right)^k \frac{1}{k!} \ & ≤ & \frac{|t|^3}{\sqrt{n}} exp\left(\frac{|t|}{\sqrt{n}}\right) →_{n → ∞} 0 \, , \end{array} \end{displaymath} which completes the proof.

Cumulative distribution functions (CDF) of $Y_5$ (black) and $Z$ a standard normal (red).

Cumulative distribution functions (CDF) of $Y₅₀$ (black) and $Z$ a standard normal (red).

Cumulative distribution functions (CDF) of $Y₅₀₀$ (black) and $Z$ a standard normal (red).

Cumulative distribution functions (CDF) of $Y₅₀₀₀$ (black) and $Z$ a standard normal (red).

CCD calibration

CCD calibration

If what I just exposed is correct, with the two (main) “noise” sources, the observations $Y$ (from a CCD pixel) follow: \[Y ∼ G \, λ + \sqrt{G^2 \, (λ+σ_R^2)} \, ε \; ,\] where $G$ is the camera gain, $σ_R^2$ is the read-out variance and $ε$ is a standard normal rv. The values of $G$ and $σ_R^2$ are specified by the manufacturer for each camera, but experience shows that manufacturers tend to be overoptimistic when it comes to their product performances—they can for instance give an underestimated $σ_R^2$. Its therefore a good idea to measure these parameters with calibration experiments. Such calibration experiments are also the occasion to check that our simple model is relevant.

• Our problem becomes: How to test $Y ∼ G \, λ + \sqrt{G^2 \, (λ+σ_R^2)} \, ε$ ? Or how to set different values for $λ$?
• Let’s consider a pixel of our CCD “looking” at a fixed volume of a fluorescein solution with a given (and stable) concentration. We have two ways of modifying λ :
  • Change the intensity $i_e$ of the light source exciting the fluorophore.
  • Change the exposure time $τ$.

We can indeed write our $λ$ as: \[λ = φ v c i_e τ \, ,\] where

• $v$ is the solution’s volume “seen” by a given pixel,
• $c$ is the fluorophore’s concentration,
• $φ$ is the quantum yield.

In practice it is easier to vary the exposure time τ and that’s what was done in the experiments described next… Question: Can you guess what these experiments are?

Sebastien Joucla and myself asked our collaborators from the Kloppenburg lab (Cologne University) to:

• choose 10 exposure times,
• for each of the 10 times, perform 100 exposures,
• for each of the 10 x 100 exposures, record the value $y_ij$ of the rv $Y_ij$ of CCD’s pixel $i,j$.

We introduce a rv $Y_ij$ for each pixel because it is very difficult (impossible) to have a uniform intensity ($i_e$) and a uniform volume ($v$) and a uniform quantum yield ($φ$). We have therefore for each pixel: \[Y_i,j ∼ G \, p_i,j τ + \sqrt{G^2 \, (p_i,j τ+σ_R^2)} \, ε_i,j\; ,\] where $p_i,j = c φ_i,j v_i,j i_e,i,j$.

• If our model is correct we should have for each pixel $i,j$, for a given exposure time, a mean value: \[\bar{y}_i,j = \frac{1}{100} ∑_k=1¹⁰⁰ y_i,j,k ≈ G \, p_i,j τ \]
• and a variance: \[S_i,j^2 = \frac{1}{99} ∑_k=1¹⁰⁰ (y_i,j,k-\bar{y}_i,j)^2 ≈ G^2 \, (p_i,j τ+σ_R^2) \; .\]
• The graph of $S_i,j^2$ vs $\bar{y}_i,j$ should be a straight line with slope $G$ ordinate at 0, $G^2 σ_R^2$.

The first exposure of 10 ms (experiment performed by Andreas Pippow, Kloppenburg Lab., Cologne University).

CCD calibration: Checking the assumptions

• The data are going to be analyzed as if the $Y_i,j,k$ were IID, but they were sequentially recorded. It is therefore strongly recommended to check that the IID hypothesis is reasonable.
• The small example of the next figure shows that there are no (obvious) trends.
• We must also check the correlation function.

Counts time evolution for three neighboring pixels (10 ms exposure time).

• If the $Y_i,j,k$ are not IID we expect a more or less linear trend—due to bleaching of the dye.
• Rather than looking at each individual pixel sequence like on the previous slide, we can fit the following linear model model to each pixel: $$Y_i,j,k = β_0 + β_1 k + σ ε_i,j$$ where the $ε_i,j \stackrel{IID}{∼} \mathcal{N}(0,1)$, and check if $β_1$ can be reasonably considered as null; while a trend due to bleaching would give a negative $β_1$.
• Without a trend, the theoretical distribution of $\hat{β}_1 / \hat{σ}_{β_1}$ —$\hat{β}_1$ is the estimate of $β_1$ and $\hat{σ}_{β_1}$ its estimated standard error—is a Student’s t distribution with 97 degrees of freedom.
• Applying this idea to the central pixel of the previous slide we get…

file:figs/middle_adu_seq_fit_fig.png

We get $\hat{β}_1 = 0.032$ and a 95 % conf. int. for it is: $[-0.018,0.082]$.

Empirical density of $\hat{β}_1 / \hat{σ}_{β_1}$ from all exposures (in black), theoretical one (t with 97 df) in orange.

• We now look for potential correlations between recording from different pixels.
• We do that by computing the empirical correlation between pixels $(i,j)$ and $(u,v)$.
• We get the empirical mean at each pixel (for a given exposure time) that is: $\overline{Y}_ij = (1/K) ∑_k=1^K Y_ijk$.
• We get the empirical variance: $S^2_ij = 1/(K-1) ∑_k=1^K (Y_ijk-\overline{Y}_ij)^2$.
• We then obtain the normalized signal or standard score: $N_ijk = (Y_ijk-\overline{Y}_ij)/\sqrt{S^2_ij}$.
• The correlation coefficient is then: $ρ(ij,uv) = 1/(K-1) ∑_k=1^K N_ijk N_uvk$.
• Under the null hypothesis, no correlation, $ρ(ij,uv) ∼ \mathcal{N}(0,1/K)$.

Empirical density in black (all exposure times using nearest neighbor pixels), theoretical one, $\mathcal{N}(0,0.01)$, in orange.

CCD calibration: again

We wrote previously :

• If our model is correct we should have for each pixel $i,j$, for a given exposure time, a mean value: \[\bar{y}_i,j = \frac{1}{100} ∑_k=1¹⁰⁰ y_i,j,k ≈ G \, p_i,j τ \]
• and a variance: \[S_i,j^2 = \frac{1}{99} ∑_k=1¹⁰⁰ (y_i,j,k-\bar{y}_i,j)^2 ≈ G^2 \, (p_i,j τ+σ_R^2) \; .\]
• The graph of $S_i,j^2$ /vs/ $\bar{y}_i,j$ should be a straight line with slope $G$ ordinate at 0, $G^2 σ_R^2$.

$S_i,j^2$ vs $\bar{y}_i,j$

We do see the expected linear relation: $\mathrm{Var}[ADU] = G^2 σ_R^2 + G \mathrm{E}[ADU]$.

Linear fit

The heteroscedasticity (inhomogeneous variance) visible on the graph is also expected since the variance of a variance for an IID sample of size $K$ from a normal distribution with mean $μ$ and variance $σ^2$ is: $$\mathrm{Var}[S^2] = \frac{2σ^4}{(K-1)} \; .$$

• This means than when we do our linear fit, $$y_k = a + b x_k + σ_k ε_k \, ,$$ we should use weights.
• Here $$x_k = \overline{ADU}_k \quad y_k = \mathrm{Var}[ADU]_k \, ,$$ $$b = G \quad a = G^2 σ_R^2 \, .$$
• It’s easy to show that the least square estimates are: $$\hat{a} = \frac{1}{Z} ∑_k \frac{y_k-\hat{b} x_k}{σ_k^2} \quad \text{where} \quad Z = ∑_k \frac{1}{σ_k^2}$$ and $$\hat{b} = \left(∑_k \frac{x_k}{σ_k^2} \left(y_k - \frac{1}{Z}∑_j \frac{y_j}{σ_j^2}\right)\right) / \left(∑_k \frac{x_k}{σ_k^2}\left(x_k - \frac{1}{Z}∑_j \frac{x_j}{σ_j^2}\right)\right) \, .$$
• We don’t know $σ_k$ but we have an estimation: $\hat{σ}_k^2 = \mathrm{Var}[S_k^2]$ we can “plug-in” this value to get our weights.

We have here $\hat{G} = 0.14$ and $\hat{σ}_R^2 = 290$.

Checking the fit

Some remarks

• When we use a linear regression, we are (implicitly) assuming that the “independent” variable, here $\overline{ADU}_k$, is exactly known.
• This was clearly not the case here since $\overline{ADU}_k$ was measured (with an error).
• We could and will therefore refine our fit.

Error propagation and variance stabilization

Error propagation

• Let us consider two random variables: $Y$ and $Z$ such that:
• $Y ≈ \mathcal{N}(μ_Y,σ^2_Y)$ or $Y ≈ μ_Y + σ_Y \, ε$
• $Z = f(Y)$, with $f$ continuous and differentiable.
• Using a first order Taylor expansion we then have:\[ \begin{array}{lcl} Z & ≈ & f(μ_Y + σ_Y \, ε) \ & ≈ & f(μ_Y) + σ_Y \, ε \, \frac{d f}{d Y}(μ_Y) \end{array}\]
• $\mathrm{E}Z ≈ f(μ_Y) = f(\mathrm{E}Y)$
• $\mathrm{Var}Z ≡ \mathrm{E}[(Z-\mathrm{E}Z)^2] ≈ σ^2_Y \, \frac{d f}{d Y}^2(μ_Y)$
• $Z ≈ f(μ_Y) + σ_Y\left| \frac{d f}{d Y}(μ_Y)\right| \, ε$

Variance stabilization: Theory

• For our CCD model we have (for a given pixel): \[Y ∼ G \, λ + \sqrt{G^2 \, (λ+σ_R^2)} \, ε = μ_Y + \sqrt{G \, μ_Y + G^2 \, σ_R^2} \, ε \, .\]
• Then if $Z = f(Y)$ we get: \[Z ≈ f(μ_Y) + \mid f’(μ_Y) \mid G \sqrt{μ_Y / G+σ_R^2} \, ε\]
• What happens then if we take: $f(x) = 2 \, \sqrt{x/G + σ_R^2}\;$?
• We have: \[f’(x) = \frac{1}{G \sqrt{ x / G + σ_R^2}}\]
• Leading to: \[Z ≈ 2 \, \sqrt{μ_Y / G + σ_R^2} + ε\]

Variance stabilization: Example

Application

Back to where we started

ADU counts (raw data) from Fura-2 excited at 340 nm. Each square corresponds to a pixel. 25.05 s of data are shown. Same scale on each sub-plot. Data recorded by Andreas Pippow (Kloppenburg Lab. Cologne University).

Quick ROI detection: Motivation

• After variance stabilization: $Z_i,j,k = 2 \, \sqrt{ADU_i,j,k / G + σ_R^2}$, the variance at each pixel $(i,j)$ at each time, $k$, should be 1.
• If a pixel contains no dynamical signal—that is nothing more than a constant background signal—the following statistics: \[RSS_i,j ≡ ∑_k=1^K (Z_i,j,k - \overline{Z}_i,j)^2 \quad \mathrm{with} \quad \overline{Z}_i,j ≡ \frac{1}{K} ∑_k=1^K Z_i,j,k\] should follow a $χ^2$ distribution with $K-1$ degrees of freedom.
• We could therefore compute the values of the complementary cumulative distribution function of the theoretical $χ_K-1^2$ distribution:\[1 - F_{χ_{K-1}^2}(RSS_i,j)\] and look for very small values—that is very small probabilities—(using a log scale helps here).

Grey scale map of $log\left(1 - F_{χ_{K-1}^2}(RSS_i,j)\right)$

Pointwise time course estimation

• We are going to be (very) conservative and keep as our ROI the pixels having an $log\left(1 - F_{χ_{K-1}^2}(RSS)\right) ≤ -300$.
• We are then left with 12 pixels.
• We are going to model the fluorescence intensity of each of these pixels by: \[S_i,j(t) = φ_i,j \, f(t) + b \; ,\] where $f(t)$ is a signal time course to all pixels of the ROI, $φ_i,j$ is a pixel specific parameter and $b$ is a background fluorescence assumed identical for each pixel.
• The time $t$ is in fact a discrete variable, $t = δ \, k$ ($δ$ = 150 ms) and we are seeking a pointwise estimation: $\{f_1,f_2,\ldots,f_K\}$ ($K$ = 168) where $f_k = f(δ \, k)$.
• We end up with 12 ($φ_i,j$) + 168 ($f_k$) + 1 ($b$) = 181 parameters for 12 x 168 = 2016 measurements.
• We need to add a constraint since with our model specification: \[S_i,j,k = φ_i,j \, f_k + b \; ,\] we can multiply all the $φ_i,j$ by 2 and divide all the $f_k$ by 2 and get the same prediction.
• We are going to set $f_k=1$ for the first 5 time points (the stimulation comes at the 11th) and our pointwise estimation relates to what is usually done with this type of data, $Δ S(t) / S_0$ (where $S_0$ is a baseline average) through: \[Δ S(t) / S_0 = \frac{S(t) - S_0}{S_0} = f(t) - 1 + \mathrm{noise}\, .\]
• Notice that no independent background measurement is used.
• With variance stabilization we end up minimizing: \[RSS\left(b,(φ_i,j),(f_k)_{k=6,\ldots,168}\right) = ∑_{(i,j) ∈ \mathrm{ROI}} ∑_k=1¹⁶⁸ \left(Z_ijk-F_ijk\right)^2 \, ,\] where \[Z_ijk = 2 \, \sqrt{ADU_ijk/\hat{G}+\hat{σ}_R^2}\] and \[F_ijk = 2 \, \sqrt{φ_i,j \, f_k + b + \hat{σ}_R^2}\, .\]
• If our model is correct we should have: \[RSS\left(\hat{b},(\hat{φ}_i,j),(\hat{f}_k)_{k=6,\ldots,168}\right) ∼ χ^2_{12 × 168 - 175}\, .\]
• The method also generates confidence intervals for the estimated parameters.

Technical details

• To solve this 175 dimensional optimization problem in a reasonable time (< 5 s) in C with the GSL library we use the Levenberg-Marquardt Algorithm performing nonlinear least-squares minimization.
• To improve numerical behavior we work with the log of the parameters, since all parameters are positive.
• Giving all the details would be at least as long as the present talk, but they are fully disclosed in the source file of this talk that can be found on Github: https://github.com/christophe-pouzat/ENP2017.

Time course estimate

Data and fit

Data and fit after variance stabilization. The RSS is 1952 giving a probability of 0.96 (a bit large).

Thanks

This work was done in collaboration with:

• Sebastien Joucla
• Romain Franconville
• Andeas Pippow
• Simon Hess
• Peter Kloppenburg

Thank you for your attention!