find peaks from response signal

Question

My subject is to model the response of micro-organisms to the pollution of the media (one specific pollutant). We analyse the response (production of a substance) through time that is linked to the dose of pollutant.

We know the kinetics of the response to a pulse/peak of pollutant, it has this shape:

The actual response that we observe (every hour) is the summation of several pulses (peaks of the pollutant through time), it is much more noisy and I want an algorithm capable of finding when the probable pulses of pollutant happened, and estimate what their intensity are.

I am thinking that it is close to some algorithms of mixture models (but not Gaussian of course!), or some Bayesian or deconvolution algorithm but can't find what method would fit my problem.

PS: I work with R.

Answer 1

Deconvolution is the best way to categorize your problem if

1. There is a standard pulse, called the point spread function (PSF), whose response over time is known
2. Response scales linearly with pulse amplitude
3. Overlapping pulse responses combine additively.

If all these assumptions hold, you have a standard linear deconvolution problem. Reliable algorithms for solution exist and used routinely in signal and image processing. Unfortunately R seems to offer much less than MATLAB for deconvolution. This implementation of the Richardson-Lucy algorithm was the best I could come up with. However, if you are analyzing count data and the response peaks are not too overlapped, R-L deconvolution could give good results.

If in addition there are only a few well-separated pulses, you can use sparse regularized deconvolution. This would likely give higher-quality, more reliable results than any other approach. There is a well-developed theory* that says when it is possible to recover the underlying pulses from noisy and limited pulse response data. It may also work if some of the assumptions above are violated. See it in action in this MATLAB demo, where you can see the recovery is very accurate.

Generally, mixture modeling algorithms assume uncertainty in the PSF and attempt to infer it from the data, a form of blind deconvolution. The results will be lower quality and less robust if you make the algorithm guess the PSF when you know what it actually is. Unfortunately I have never seen a mixture modeling package in R that allows you to help the algorithm by specifying the known PSF. However if there is uncertainty in your PSF, you may want to consider the mixtools R package.

*From the field of sparse modeling and compressive sensing

Answer 2

If you know the formula of the response to a single pulse/peak, you could do a deconvolution of the observed response. You could use the approach described in this article http://www.pnas.org/content/109/22/8456.long. Very different field, but similar problem I think (look at their fig. 2); they also provide the code in R to run the analysis as supplemental material.

Note that you will have to specify in advance the possible locations in time of the pulses (e.g., assuming that a pulse might occur once every 20 minutes for example) and therefore their number. The optimization algorithm then will give you the intensity of these pulses. You could try with different spacing in time of the pulses locations, but be careful that if you have many pulses (very close in time one to another) the optimization might become unstable and highly dependent on the initial values (you should repeat the process many times, with different initial values, and check that you get stable estimates of pulse intensities).

Answer 3

I don't know an algorithm allowing to do that.

However, your problem seems to be related to the control theory. From my point of view, I will model the distribution of the polluant through time using a two compartment model. For instance, you can use the ODE model:

$\dot{R} (t) = − ω \times (R (t) − X (t))$

$\dot{X} (t) = − ω \times (X (t) − P (t))$

where $P$ is your pulse ($P(t)=P$ if $t=t_p$ and $0$ otherwise) and $R$ your response (polluant concentration). From this, it can be shown that the peak of concentration is reached at

$t_{max}=1/ω+t_p$

where $t_p$ is the time of the pulse. (see (p.113): https://tel.archives-ouvertes.fr/tel-00968903/document)

As a result, you can linked the time of the peak to the time of the pulse. For instance, following your graph, if the peak is reached after three hours, then $1/ω=180$ if you use minutes. Next, you can try to find the optimal value of $P$ minimizing the errors between the ODE model and your data...