 Research
 Open Access
 Published:
Novel EM based ML Kalman estimation framework for superresolution of stochastic threestates microtubule signal
BMC Systems Biology volume 12, Article number: 112 (2018)
Abstract
Background
Recent research has found that abnormal functioning of Microtubules (MTs) could be linked to fatal diseases such as Alzheimer’s. Hence, there is an imminent need to understand the implications of MTs for disease diagnosis. However, studies of cellular processes like MTs are often constrained by physical limitations of their data acquisition systems such as optical microscopes and are vulnerable to either destruction of the specimen or the probe. In addition, study of MTs is challenged with nonuniform sampling of the MT dynamic instability phenomenon relative to its timelapse observation of the cellular processes. Thus, the above caveats limit the overall period of time that the MT data can be collected, thereby causing limited data availability scenario.
Results
In this work, two novel superresolution frameworks based on Expectation Maximization (EM) based Maximum Likelihood (ML) estimation using Kalman filters (MLK) technique are proposed to address the issues of nonuniform sampling and limited data availability of MT signals. The proposed MLK methods optimizes prediction of missing observations in the MT signal through information extraction using correlationbased patch processing and principal component analysis based mutual information. Experimental results prove that the proposed MLKbased superresolution methods outperformed nonlinear interpolation and compressed sensing methods.
Conclusions
This work aims to address limited data availability and data/observation loss incurred due to nonuniform sampling of biological signals such as MTs. For this purpose, statistical modelling of stochastic MT signals using EM based ML driven Kalman estimation (MLK) is considered as a fundamental framework for prediction of missing MT observations. It was experimentally validated that the proposed superresolution methods provided superior overall performance, better MT signal estimation using fewer samples, high SNR, low errors, and better MT parameter estimation than other methods.
Background
Research on Microtubules (MTs) have recently garnered a great level of interest due to its anomalous functioning being associated with the onset of several lethal diseases including Alzheimer’s disease [1], Parkinson’s disease [2], and various forms of cancer [3]. Essentially, MTs are intracellular polymers made of tubulin protein dimers [4] which are found in all eukaryotic cells, and play major roles in many intracellular activities such as trafficking, mitosis, cell motility [4], and chromosome segregation [5]. Under fixed external conditions, it was first noted that the MTs reached a steady state while randomly switching between polymerization (growth) and depolymerization (shrinking) states [4]. This unique phenomena of MTs randomly switching between the two states was named as “dynamic instability” by Mitchison and Kirschner [4, 6–8]. It was later reported in [9, 10], that the MTs were in fact switching spontaneously between three states, namely, growth (g), pause (p) and shrinkage (s). This threestate dynamic instability model involves eight parameters: six transition rates between the states of growth, pause and shrinkage and growth and shrinkage velocities v_{g} and v_{s}. The average length and mean lifetime of a MT depends on the threshold quantity of these parameters.
Optical microscopes have been utilized for decades as a mainstay for data collection of various cellular processes, such as MTs. Even after modern technological advances in the area of microscopy like single molecule sensitivity, and frame rates in microseconds, optical microscopes are still vulnerable in either damaging the specimen of interest, or the probe during the course of experiment. Moreover, the probe can only be illuminated for a certain period of time, making data collection an arduous task, and data availability is renderedsparse with respect to the time scale of intracellular processes [11, 12]. Especially, in case of MTs, where its high temporal resolution is imperative to estimate the dynamic instability parameters, understand the MT behavior, and other cellular activities. In particular, estimation of MT dynamic instability parameters is crucial to understand the state of the MT system and study its relationships to the occurrence of fatal diseases like Alzheimer’s. Consequently, it signifies the need for stochastic methods that can analytically overcome the challenges such as scarce data availability, and equipment limitations through better signal reconstruction and estimation techniques.
Many statistical methods for modeling stochastic signals have been suggested in literature such as hidden Markov models (HMMs) for capturing the underlying signal variability [13]. As a variation of traditional HMMs, many waveletbased HMMs have evolved because of their capability to model realworld nonGaussian signals [14, 15]. Compressed sensing (CS) [16, 17] and other subspace learning methods [18] have also been proposed to reconstruct signals using fewer samples than the Nyquist rate. In literature, it has also been proposed to use maximum likelihood (ML) technique to estimate parameters of a linear dynamic system from the observed data [19, 20]. Typically, this involves formulation of a timevarying Kalman predictor with a likelihood function to minimize the prediction error. In this case, a convergence assumption is made for sufficiently large observations which transforms the minimization of likelihood function to a nonlinear programming problem [21–23]. As an alternative, it was proposed to use ExpectationMaximization (EM) based maximum likelihood (ML) estimation in Kalman filters with the assumption that the state is now observable. This assumption simplifies the convergence criteria by eliminating the need for a large observation limit and facilitates to solve cases with missing observations [24]. Kalman filters have been proved to be optimal in minimizing mean square error sense, under the assumption that all noise is Gaussian. Therefore, the use of EM based ML estimation in Kalman filters as a means to study MTs is very relevant, since we want to minimize the residual error between the original MT signal, and the predicted MT signal.
In this paper, the MT signal is modelled as a threestate stochastic random evolution signal on which nonuniform sampling is performed to emulate the data loss in realworld scenario [11, 12]. Our motivation is to improve prediction of the missing timelapse intracellular observations analytically; in an effort to minimize the effect of external dependencies such as spatiotemporal resolution, and experimental equipment precision. For this, we propose two novel methods for superresolution of MT signals based on nonuniform sampling, using EM based ML estimation Kalman filter for better prediction of the interpolated MT signal, which is followed by correlation coefficient (R) based patch processing (MLKR) [25] and principal component analysis (PCA)based mutual information (MI) criterion (MLKMI) for information extraction to further optimize our final predicted MT signal. We estimate the dynamic instability parameters for three states in MTs using waveletbased peak detection and compare it with our previous work using CS [16]. Our proposed methods MLKR and MLKMI gave superior overall performance compared to all other methods, better MT signal recovery and gave high SNR with low errors validating the efficacy of our approaches.
Methods
Nonuniform sampling and interpolation
Consider a stochastic input MT signal x(n) with N samples. Nonuniform sampling is performed on this MT signal x(n) to generate its downsampled low resolution version x_{l}(n). The nonuniform sampling case is specifically considered to emulate the loss of information occurring in the realworld scenario due to hardware limitations of the data acquisition system [11, 12]. Therefore, the assumption is that the only signal available at hand for data analysis is the low resolution MT signal x_{l}(n). To accommodate the missing values due to nonuniform sampling, a nonlinear interpolation scheme such as piecewise linear interpolation is then used on the low resolution MT signal x_{l}(n) to get its corresponding high resolution interpolated version x_{h}(n). From this high resolution interpolated frame we can generate a dictionary of d new low resolution frames through interlaced nonuniform sampling, where each frame is derived with a constant shift factor ε_{i} from the first frame as below:
where i=2,…,d+1, and x_{l1}(n)=x_{l}(n), mod=modulo operator. For simplicity, we choose ε_{i}=i+1, and d is the resolution factor and it also denotes the level of downsampling performed. For e.g. d=2 implies that the signal was downsampled by a factor of 2. s is the number of samples in each frame, calculated as \( s= \left \lfloor \frac {N}{d}\right \rfloor \). For consistency, nonlinear interpolation using piecewise linear interpolation method is performed on the d new low resolution MT signals x_{li}(n) to get their corresponding high resolution MT signals y_{hi}(n).
Expectationmaximization based maximum likelihood estimation of a stochastic MT system
From the previous section, it is assumed that we are only given the low resolution MT signal and that nonlinear interpolation is used as a basis to provide a first estimate for missing signal values and yield the corresponding high resolution MT signals. Thus, the problem of estimation of missing observations of a stochastic MT signal from few known measurements can be formulated as a d series of Maximum Likelihood (ML) parameter estimation for a timevarying linear dynamic MT system from observations at hand y_{hi}(n). Our goal is to achieve the best approximation \(\hat {\boldsymbol{x}}_{hi}(n)\) of the unknown input signal (or state vector) x(n) from observed y_{hi}(n). Let the set of observations Y={y_{h1}(n),y_{h2}(n),…,y_{hd+1}(n)} be generated by a linear dynamic MT system. Then the timevarying Kalman filter equations describing the system can be formulated as below (using notation x_{n}=x(n)):
The F matrix represents the state at the previous time step n to the state at the current step n+1 in the absence of process noise. The matrix H gives the relationship between the state x_{n} to the measurement y_{n}. Both process noise w_{n} and measurement noise v_{n} are assumed as uncorrelated zero mean white Gaussian noise vectors with covariances:
Where T denotes transpose of a matrix, Q represents process noise covariance and R for measurement noise covariance. In practice, the matrices F, H, Q and R are subject to change with each time step or measurement (See [22] for detailed information). It is also assumed that the initial conditions for state x_{0} is Gaussian with a known mean and covariance (μ_{0},Σ_{0}). Theerefore, the ML estimates of the unknown parameters θ in F, H, Q, R, can be obtained by minimizing the negative log likelihood or as in [19, 20]:
Where e_{n}, \(\Sigma _{e_{n}}\) are the prediction error and its covariance, and can be obtained from Kalman filter time and measurement update equations as in [20, 22] below:
Time update:
Measurement update:
where K is the Kalman gain. The function of the measurement update is to adjust the projected estimate by the actual measurement at that time n, whereas time update serves the purpose of prediction of the current state estimate \(\hat { \boldsymbol {x}}_{n+1n}\) ahead of time (indicated through future time instance n+1). This recursive update is performed with the objective of minimizing (4) with respect to unknown parameter θ to estimate the missing observations. In this work, we have used Expectation Maximization (EM) algorithm to find the ML estimates as in [22–24], by the assuming that the state is now observable and can be denoted as X={x_{h}(0),…,x_{h}(n)}, hence the problem in (4) can be reformulated as:
Using fixed interval form of Kalman filter (RTS smoother), the new system recursion updates can be computed as in [22]:
Forward Recursions:
Backward Recursions:
Note that the formulation of d series estimation of the state vector x_{n} is made with the assumption that available data is the observed vector y_{n} and it provides d unique estimated MT signals \(\hat { \boldsymbol {x}}_{hi}(n)\). In each iteration, EM algorithm computes the datasufficient statistics in recursions (8), (9) and estimates previous model parameters (Estep). New system parameters are obtained from these statistics in the maximization step (Mstep).
Correlation coefficientbasedpatch processing
Typically, employing correlationbased approaches are a common practice in superresolution for images. This is because correlation coefficient (R) is highly effective in detection and extraction of identical information present in the low resolution images. Hence, inspired by the proficiency of correlation coefficient in information extraction, we have implemented a novel correlation coefficient (R) based patch processing technique to find an identical match between similar signal patches (segments) present in the d+1 Kalmanpredicted MT signals \(\hat { \boldsymbol {x}}_{hp}(n)\) and achieve further optimization to yield an enhanced final prediction for MT signal. This is accomplished by using a sliding window w=4(w=1×4) across the pairwise MT signals. The pairwise R comparisons ensures that the inconsistencies such as artificial low/high frequencies introduced during signal sampling or processing are removed, and maximum information content is extracted from the aliased signals. In order to define a patch as similar or identical the R threshold parameter is chosen as R>0.90. The MT signal patches that satisfy this condition as in (10) are then included in the final predicted signal \(\hat {\boldsymbol {x}}_{f}(n)\). An error minimization condition is also imposed such that the selected patch must not only satisfy R>0.90, but also minimize the error between the signals as below:
where p, q=1,2,…,d. For d Kalmanpredicted MT signals \(\hat { \boldsymbol {x}}_{hp}(n)\), we will need ^{d}C_{2} signal comparisons to be made. All missing values in the final MT signal \( \hat {\boldsymbol {x}}_{f}(n)\) (for patches < threshold R) are computed by taking the mean of the d signals values available for that particular time instance. This step is done as a tradeoff to minimize the prediction error for missing values in the original low resolution MT signal. The final MLK−Rpatch processed MT signal \( \hat {\boldsymbol {x}}_{f}(n)\) had better SNR and lower errors than Kalmanpredicted (MLK) (\(\hat { \boldsymbol {x}}_{hi}(n)\)) and interpolated (NLI) (y_{hi}(n)) signals. This final reconstructed MT signal \(\hat {\boldsymbol {x}}_{f}(n)\) was used for estimation of the threestate MT dynamic instability parameters in the wavelet domain using peak detection as in our prior works [16, 26].
Principal component analysisbased mutual information criterion
Principal component analysis (PCA) is a widely used unsupervised learning technique for information extraction and dimensionality reduction of data in multidisciplinary applications. The principle behind PCA is that it captures the directions of maximum variance (information) in the data, where these directions of maximum variance represent the principal components or eigenvectors of the data [27]. The significance of PCA includes: 1) the principal components are uncorrelated, thus they provide elimination of redundancy in the data. 2) Most of the important information is usually compacted in the first few eigenvectors, thereby providing data compaction and effective dimensionality reduction of data. In this work, we apply PCA to extract information present in the individual d+1 MLK predicted MT signals. Given a data x, it can be described as a linear combination of its eigenvectors using PCA as:
where A denotes the new principal component scores, B gives weight vectors such that the linear combination of xB maximizes the variance contained in A and data x can be recovered as:
Mutual information (MI) is a wellknown information extraction criterion that measures the similarity of information content or mutual dependence between two random variables. It was first conceptualized by Shannon [28], since then it been widely adopted in a variety of applications, especially in the area of biology [29, 30]. The mutual information I between two discrete random variables X and Y can be described as in [29]:
where p_{X,Y}(x,y) stands for the joint probability density function of X and Y, and p_{X}(x) and p_{Y}(y) denote the marginal probability density functions of X and Y respectively. Intuitively, if random variables X and Y are independent, then their MI I(X,Y)=0, whereas if they have perfect dependence then their MI tends to infinity i.e., I(X,Y)→∞. As defined by Shannon in [28], mutual information I(X,Y) can also be defined in terms of entropy as [28]:
Where H(X) denotes entropy of random variable X, H(Y) represents entropy of random variable Y and H(X,Y) describes the joint entropy of random variables X and Y. The above definition for MI can be interpreted in terms of entropy and probability density functions of the corresponding random variables by substituting Eqs. (13) in (16) to yield:
Thus, in this work, we apply MI to determine the information dependence between the principal components of the d+1 predicted MT signals. Greater MI signifies more common information content between the signals. Since we are trying to refine our estimation of the missing values in the MT signal, it is pertinent to retain the common information in the MT predicted signal. Lesser MI could imply a bad signal prediction and therefore should be excluded from analysis to improve final MT signal prediction \( \hat {\boldsymbol {x}}_{f}(n)\). The first two unique principal components p corresponding to the MI pairs that give us the maximum MI can be found as below:
where i,j∈{1,d+1} and I(x_{ai},x_{aj}) denotes the MI between pairs of principal components ai and aj of data x. The average of p principal components selected from our PCAMI criterion (MLKMI) is then computed to reconstruct our final MT signal \( \hat {\boldsymbol {x}}_{f}(n)\) using (12). Figure 1 illustrates the block diagram of the proposed superresolution framework for MT signal prediction of missing observations.
Wavelet transforms
Wavelet transforms are commonly used in a wide variety of biomedical applications like for compression of EEG and ECG signals. Wavelets are specially preferred for modelling stochastic signals due to its attractive benefits such as sparsity, compression, inherent denoising and simultaneous timefrequency resolution of signals. Given a 1Ddiscretetime signal \( \hat {\boldsymbol {x}}_{f}(n) \), it can be represented using translations of basis mother wavelet function ψ(n) and a lowpass scaling function ϕ(n) in the time and frequency domain as in ([15, 31]) by:
where \( \hat {\boldsymbol {x}}_{f}(n) \in L_{2} (R)\), k,j∈Z. k stands for time (translation) and j represents the frequency (scale) respectively [31] and the basis functions are:
where the first term in (17) corresponds to the coarse resolution while the second term represents the detail (or wavelet) resolution of the signal. \(c_{j_{0}}(k)\) and d_{j}(k) are the corresponding approximation (or scale) and detail (or wavelet) coefficients at scale j, respectively. These coefficients can be calculated as below [31]:
Note that j_{0} is an arbitrary starting scale. In this paper, the maximum scale j also represents the wavelet decomposition levels. We later perform peak detection of MTs in the wavelet domain using the energy packing density (EPE) criterion [32] to estimate the threestate MT parameters. In particular, simultaneous timefrequency resolution is the key wavelet property that will be used for detection of the threetransition states in MTs.
Trichotomous Markov Noisebased threestates random evolution model
Consider the final predicted MT signal \(\hat {\boldsymbol {x}}_{f}(n)\) that is undergoing the dynamic instability phenomena by randomly switching between the three transition states of growth (g), pause (p) and shrinkage (s). We formulate a threestates random evolution model for the study of dynamic instability process in MTs using Trichotomous Markov Noise (TMN), which is a threestate stochastic process. Hence, the threestates random evolution model for the dynamic instability of MTs has eight statesspace parameters, including six transition rates between states of growth, pause and shrinkage denoted by f_{ij}, where i,j∈{g,p,s} and two states of growth and shrinkage rates, represented by v_{g} and v_{s} [16]. The mean length and lifetime of a MT in the threestates random evolution model depends on the threshold quantities f_{gs}, f_{sg}, v_{g}, and v_{s} in the formula for V in (20) [33]. If the quantity, V is positive, the MTs tend to shrink more than they grow, and the MTs will have a finite mean length and mean lifetime. Otherwise, on average, they tend to grow forever. The derivation of equations for threestates random evolution model for MTs are provided in [16]:
Where V gives the equilibrium point of the system and L denotes the average length of the MT signal.
Results and discussion
In this section, we experimentally validate the effectiveness of our proposed methods for superresolution of MT signals. The MT data used in this work are results of experiments on MTs (composed of purified αβ_{II} isotopes from bovine brain tubulin) performed by O. Azarenko, L. Wilson and M.A Jordan at the University of California, Santa Barbara. Tubulin proteins were first purified from the bovine brain and then seeded to polymerize at 37°C. The growth and shrinkage dynamics of individual purified MTs were then recorded at their plus ends using the differential interference contrast video microscopy. Data points representing MT lengths were collected at 2−6s intervals. MT lengths were analyzed using the Real Time Measurement program. Growth and shrinkage rates were calculated by leastsquares regression analysis of the data points. Growth and shrinkage thresholds are set to an increase in length by 0.2μm at a rate of 0.15μm/min and a decrease in length by 0.2μm at a rate of 0.3μm/min, respectively. Any length changes equal to or less than 0.2μm over the duration of six data points were considered attenuation phases (phases in which length changes were below the resolution of the microscope). It should be noted that the experimental detection limit for length changes corresponds to about 400−800 tubulin dimers. The supplied data, however, was in the form of a hard copy graphs, that they were scanned and then digitized using the software “DigitizeIt” (http://www.digitizeit.de/) [34].
The MT data used in this study is ABII (αβ_{II}) data and the number of samples present in this original MT signal was N=165. In this study, nonuniform sampling was performed on the original MT signal by a factor of d=3 to emulate the inadequate MT data collection limitation in physical systems. Figure 2 shows the original MT signal. The red points indicate the nonuniformly chosen samples, while blue points indicate the samples chosen through uniform sampling by a factor d=3. Note that the blue points were more evenly distributed, unlike red points that tend to be clustered around certain areas indicating a nonuniform distribution/sampling across the MT signal.
In this paper, we introduce two novel methods for enhanced MT signal prediction in the new superresolution framework for MTs. Both the proposed methods have nonuniform sampling performed on the MT signals to emulate realworld data loss scenario. This is followed by EMbased ML Kalman prediction (MLK) of the MT signals as a basis framework for superresolution. The first proposed enhanced prediction method involves using an image processinginspired correlation (R)basedpatch processing (MLKR) method to further optimize the final MT predicted signal by extraction of information through identification of similar patches from the dictionary of MLKpredicted MT signals. Second method involves PCA based MI criterion (MLKMI) method that performs extraction of similar information present in principle components of the MLKpredicted MT signals using MI criterion for elimination of data redundancy and further enhancement of the predicted MT signal. Comparison analysis of the proposed superresolution methods is performed with respect to two other methods, namely, nonuniform sampling of MT signal followed by nonlinear (piecewise linear) interpolation technique (NLI) and our previous work using compressed sensing (CS) to reconstruct MT signals with fewer samples than the Nyquist rate [16]. Since, the nonuniformly sampled low resolution MT signal was sampled at a factor of d=3, this implies that it has \(s= \left \lfloor \frac {N}{d}\right \rfloor = \left \lfloor \frac {165}{3}\right \rfloor =55\) samples. This is equivalent to CS subrate of \( \frac {s}{N}= \frac {55}{165}= 0.33\). Therefore, we included results of CS rate =0.3 (CS0.3) from our prior work [16] for comparison.
As discussed under Methods section, after nonuniform sampling is performed to yield the resultant low resolution MT signal, its corresponding high resolution MT signal is then generated using nonlinear interpolation techniques like piecewise interpolation in our case. An interlaced nonuniform sampling of this resultant high resolution interpolated MT signal is carried out to generate d low resolution MT signals and their corresponding high resolution MT signal versions y_{hi}(n) this is denoted as the NLI method. Subsequently, EM based ML estimation using Kalman filters is applied to these high resolution MT signals to achieve a better prediction of the missing MT values (MLK method). All the methods were followed by peak detection in wavelet domain to estimate the threestate dynamic instability parameters for MTs. Figures 3 and 4 illustrates the final predicted MT signals across all methods, where Fig. 4 shows the average of d+1 predicted MT signals for both NLI and MLK methods. From Figs. 34, it can be inferred that our proposed methods MLKR and MLKMI had the best final predicted MT signal across all the methods. To quantitatively substantiate the effectiveness of our approach, we compute the signaltonoiseratio (SNR) and root mean square error (RMSE) statistics of the predicted MT signals from all methods in Table 1. Again, our proposed superresolution methods MLKR and MLKMI had the highest SNR and lowest RMSE of all other methods. In particular, MLKR outperformed all other methods both qualitatively and quantitatively.
The main difference between our proposed superresolution framework in this study and CS method from previous work [16] is that in the former case we assume that only the low resolution MT signal is available to us and using statistical modelling we try to predict and estimate the missing values in the original MT signal with minimum error possible. Whereas in the latter case, we assume that we have already acquired the original MT signal and want to compress it at the sensor side for storage or transmission purposes such that at the receiver side both the compressed signal and the transformation basis is known. The main reason for wide acceptance of CS is that we can preserve the signal details using a random measurement matrix as Gaussian, thereby eliminating the need for complex processing and transmission of transformation basis. But CS makes no effort in learning the underlying signal characteristics which can be critical in understanding biological processes as in case of MTs. On the other hand, statistical modeling methods like EM based ML estimation based Kalman prediction learns the signal structure and hence provides better signal prediction. This is demonstrated through Fig. 3 and Table 1, wherein MLKR and MLKMI methods had higher SNR, lower RMSE and performed better than the CS methods, due to the datalearning and information extraction processes involved. Thus, both our current work and previous work [16] makes an effort to analytically address various facets of the data acquisition problem in MTs and biomedical signals in general, such as low sample availability scenario, spatiotemporal resolution and compression of biomedical signals.
In particular, statistical modelling of MTs as a stochastic signal with the problem of recovering the original signal through fewer samples, while assuming all noise to be Gaussian makes Kalman prediction an obvious choice given its optimality in the mean squared error sense. Kalman filtering is also ideal for realtime signal processing, where the future samples can be predicted in realtime from the adaptive feedback of statistics computed from the current samples. A good prediction of MT signal is quintessential for better estimation of dynamic instability parameters of threestate MTs, understand the MT behavior and state of the MT system. Thus, our MLKR and MLKMI method has shown to provide better prediction and estimation of the unknown parameters in the MT linear dynamic system. In addition, wavelet domain MT signal processing is done to exploit benefits such as sparsity, inherent denoising and most importantly simultaneous timefrequency resolution which is very pertinent for good peak/ MT transition state detection. For consistency with our prior work [16], we have used Daubechies D2 (db2) as mother wavelet with 8 levels of signal decomposition; it also performs experimentally better than other mother wavelet families for our work. Wavelet transform is performed on the final predicted MT signal to get its corresponding sparse wavelet coefficients. Wavelet domain peak detection is then employed to detect the peaks indicating the transition points in the MT signal between the growth (g), pause (p) and shrinkage (s) states. This transition information is used to determine when and where the MTs are switching between the three states. This is where the simultaneous timefrequency resolution of wavelets plays a crucial role. Since, we want to detect switching instants (or time), implies we are looking for narrower time bins for better time resolution, which corresponds to the lowest level significant wavelet coefficients (swc). As in our prior work [16, 26], energy packing efficiency (EPE) [32] is used to determine the number of significant wavelet coefficients/peaks that need to be chosen. The EPE steps for peak detection of MTs are as follows:

Step 1: Sort the lowest level wavelet coefficients in descending order (i.e. larger significant wavelet coefficients are most likely peaks).

Step 2: Compute the total energy present in the wavelet coefficients (swc) in the lowest level. Where the total energy is calculated by:
$$ E_{TOT} = \sum{\mathbf{swc}^{2}} $$(21) 
Step 3: Fix a desired threshold to indicate the percentage of significant wavelet coefficients to be retained. In our case, we retain values such that 85% of total energy of the coefficients in the lowest level is preserved. That is:
$$ E_{TH} \geq 0.85*E_{TOT} $$(22) 
Step 4: Compute the total energy of the sorted wavelet coefficients (E_{TH}) as in (21), until the threshold condition in (22) is met.
The value of EPE threshold parameter plays a vital in detection of peaks and threetransition states in MTs. Higher the threshold parameter, lower the number of significant wavelet coefficients selected, thereby fewer peaks/ MT transition states are detected (implying loss of information, if threshold parameter is too high). Lower the threshold, higher the number of false peaks detected, thus causing an overprediction of the MT transition states. Therefore, the number of significant wavelet coefficients retained is sensitive to threshold chosen, which in turn dictates the number of potential peaks or MT transition states detected. Figure 5 shows the peaks detected in predicted MT signal using wavelet domain EPE method. The information derived from peak detection is then used to encode our final predicted signal to its nearest peak to obtain a TMN modelled threestate MT signal. The resultant encoded threestate MT signal is used to compute MT parameters like transition frequencies, transition velocities and average MT lengths as tabulated in Tables 2 and 3. The error estimates of MT parameters for all methods are given in Table 4. From Tables 2 and 4, we substantiate that overall MLKR and MLKMI had the best performance compared to other methods. Specially, our proposed method MLKR demonstrated superior overall performance, higher SNR, lower RMSE, better MT signal prediction and parameter estimation with lowest errors from fewer samples compared to interpolation and CS methods.
Conclusion
In this paper, we propose two novel frameworks for superresolution of MT signals to address the limited data availability scenario and emulate the data loss due to nonuniform sampling of biological/natural signals that we often encounter in the realworld. This work exploits the stochastic nature of MT signals through statistical modelling using EM based ML driven Kalman estimation (MLK) for better prediction of the nonuniformly sampled MT signal as a basic superresolution framework. The MLKpredicted MT signals were further optimized through information extraction using correlation (R)basedpatch processing (MLKR) and PCAbased mutual information criterion (MLKMI) methods. We perform comparison analysis of our proposed methods MLKR and MLKMI with respect to nonlinear interpolation (NLI) and CS methods. It was experimentally found that both the proposed methods MLKR and MLKMI achieved the best overall performance for MT signal estimation. Specifically, MLKR outperformed all the methods, and had better reconstruction performance using fewer samples, gave high SNR, low errors, as well as better MT parameter estimation than other compared methods. This work aims to demonstrate the effectiveness and significance of statistical modelling and data learning in MTs, and biomedical paradigm. Our goal is to provide an analytical solution to overcome the equipment/hardware fallacies that might occur during the signal acquisition process.
References
 1
Ballatore C, Lee VMY, Trojanowski JQ. Taumediated neurodegeneration in alzheimer’s disease and related disorders. Nat Rev Neurosci. 2007; 8(9):663–72.
 2
Farrer MJ. Genetics of parkinson disease: paradigm shifts and future prospects. Nat Rev Genet. 2006; 7(4):306–18.
 3
Kops GJ, Weaver BA, Cleveland DW. On the road to cancer: aneuploidy and the mitotic checkpoint. Nat Rev Cancer. 2005; 5(10):773–85.
 4
Mitchison T, Kirschner M, et al.Dynamic instability of microtubule growth. Nature. 1984; 312(5991):237–42.
 5
Barton NR, Goldstein L. Going mobile: microtubule motors and chromosome segregation. Proc Natl Acad Sci. 1996; 93(5):1735–42.
 6
Burbank KS, Mitchison TJ. Microtubule dynamic instability. Curr Biol. 2006; 16(14):516–7.
 7
Yarahmadian S, Barker B, Zumbrun K, Shaw SL. Existence and stability of steady states of a reaction convection diffusion equation modeling microtubule formation. J Math Biol. 2011; 63(3):459–92.
 8
Yarahmadian S, Yari M. Phase transition analysis of the dynamic instability of microtubules. Nonlinearity. 2014; 27(9):2165.
 9
Shaw SL, Kamyar R, Ehrhardt DW. Sustained microtubule treadmilling in arabidopsis cortical arrays. Science. 2003; 300(5626):1715–8.
 10
Ambrose JC, Wasteneys GO. Clasp modulates microtubulecortex interaction during selforganization of acentrosomal microtubules. Mol Biol Cell. 2008; 19(11):4730–7.
 11
Maciejewski MW, Stern AS, King GF, Hoch JC. Nonuniform sampling in biomolecular nmr. In: Modern Magnetic Resonance. Springer: 2008. p. 1305–11.
 12
Hyberts SG, Arthanari H, Wagner G. Applications of nonuniform sampling and processing. In: Novel Sampling Approaches in Higher Dimensional NMR. Springer: 2011. p. 125–48.
 13
Baum LE, Petrie T, Soules G, Weiss N. A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains. Ann Math Stat. 1970; 41(1):164–71.
 14
Crouse MS, Nowak RD, Mhirsi K, Baraniuk RG. Waveletdomain hidden markov models for signal detection and classification. In: Advanced Signal Processing: Algorithms, Architectures, and Implementations VII Vol. 3162. International Society for Optics and Photonics: 1997. p. 36–48.
 15
Crouse MS, Nowak RD, Baraniuk RG. Waveletbased statistical signal processing using hidden markov models. IEEE Trans Signal Proc. 1998; 46(4):886–902.
 16
Yarahmadian S, Menon V, Rezania V. On using compressed sensing and peak detection method for the dynamic instability parameters estimation for microtubules modeled in three states. In: Bioinformatics and Biomedicine (BIBM), 2015 IEEE International Conference On. IEEE: 2015. p. 417–20.
 17
Mishali M, Eldar YC. Blind multiband signal reconstruction: Compressed sensing for analog signals. IEEE Trans Signal Process. 2009; 57(3):993–1009.
 18
Vandewalle P, Sbaiz L, Vandewalle J, Vetterli M. Superresolution from unregistered and totally aliased signals using subspace methods. IEEE Trans Signal Process. 2007; 55(7):3687–703.
 19
Kalman RE, et al. A new approach to linear filtering and prediction problems. J Basic Eng. 1960; 82(1):35–45.
 20
Kashyap R. Maximum likelihood identification of stochastic linear systems. IEEE Trans Autom Control. 1970; 15(1):25–34.
 21
Shumway RH, Stoffer DS. An approach to time series smoothing and forecasting using the em algorithm. J Time Ser Anal. 1982; 3(4):253–64.
 22
Digalakis V, Rohlicek JR, Ostendorf M. Ml estimation of a stochastic linear system with the em algorithm and its application to speech recognition. IEEE Trans Speech Audio Process. 1993; 1(4):431–42.
 23
Ghahramani Z, Hinton GE. Parameter estimation for linear dynamical systems. Technical report, Technical Report CRGTR962, University of Totronto, Dept. of Computer Science. 1996.
 24
Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. Journal of the royal statistical society. Series B (methodological). J R Stat Soc Ser B Methodol. 1977;:1–38.
 25
Menon V, Yarahmadian S, Rezania V. Superresolution and em based ml kalman estimation of the stochastic microtubule signal modeled as three states random evolution. In: 2017 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). IEEE: 2017. p. 686–93.
 26
Mahrooghy M, Yarahmadian S, Menon V, Rezania V, Tuszynski JA. The use of compressive sensing and peak detection in the reconstruction of microtubules length time series in the process of dynamic instability. Comput Biol Med. 2015; 65:25–33.
 27
Jolliffe I. Principal Component Analysis. In: In International encyclopedia of statistical science. Berlin Heidelberg: Springer: 2002. p. 1094–6.
 28
Shannon CE. A mathematical theory of communication. ACM SIGMOBILE Mob Comput Commun Rev. 2001; 5(1):3–55.
 29
Pham TH, Ho TB, Nguyen QD, Tran DH, et al. Multivariate mutual information measures for discovering biological networks. In: Computing and Communication Technologies, Research, Innovation, and Vision for the Future (RIVF), 2012 IEEE RIVF International Conference On. IEEE: 2012. p. 1–6.
 30
Malladi R, Johnson DH, Kalamangalam GP, Tandon N, Aazhang B. Mutual Information in Frequency and Its Application to Measure Cross Frequency Coupling in Epilepsy. IEEE Transactions on Signal Processing. arXiv preprint arXiv:1711.01629. 2017; 66(11):3008–23.
 31
Burrus CS, Gopinath RA, Guo H, Odegard JE, Selesnick IW. Introduction to Wavelets and Wavelet Transforms: a Primer, vol. 1: Prentice hall New Jersey; 1998.
 32
Yip P, Rao K. Energy packing efficiency for the generalized discrete transforms. IEEE Trans Commun. 1978; 26(8):1257–62.
 33
Dogterom M, Leibler S. Physical aspects of the growth and regulation of microtubule structures. Phys Rev Lett. 1993; 70(9):1347–50.
 34
Rezania V, Azarenko O, Jordan MA, Bolterauer H, Ludueña RF, Huzil JT, Tuszynski JA. Microtubule assembly of isotypically purified tubulin and its mixtures. Biophys J. 2008; 95(4):1993–2008.
Acknowledgements
Nothing to declare.
Funding
Publication costs were funded through faculty startup funds of author VM offered by University of Alabama in Huntsville.
Availability of data and materials
Statistical and computational models used are fully detailed in the main text. Data will be made available upon personal requests to authors.
About this supplement
This article has been published as part of BMC Systems Biology Volume 12 Supplement 6, 2018: Selected articles from the IEEE BIBM International Conference on Bioinformatics & Biomedicine (BIBM) 2017: systems biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume12supplement6.
Author information
Affiliations
Contributions
VM and SY conceptualized this idea. VM implemented the MT statistical modelling and pertinent signal processing concepts. SY formulated MT mathematical modelling. VM, SY, and VR analyzed the results. VM and SY drafted the manuscript. All authors read and approved the final manuscript.
Corresponding author
Correspondence to Vineetha Menon.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Menon, V., Yarahmadian, S. & Rezania, V. Novel EM based ML Kalman estimation framework for superresolution of stochastic threestates microtubule signal. BMC Syst Biol 12, 112 (2018) doi:10.1186/s1291801806315
Published:
Keywords
 Superresolution
 Kalman filtering
 Expectation Maximization
 Wavelets
 Principal component analysis
 Mutual information
 Missing data