{'fraction_non_alphanumeric': 0.04132457580733443, 'fraction_numerical': 0.01751505199781062, 'mean_word_length': 4.231784918867324, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 14, 'https://': 4, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 72, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FMM: An R Package for Modeling Rhythmic Patterns in Oscillatory Systems

This paper is dedicated to the R package FMM which implements a novel approach to describe rhythmic patterns in oscillatory signals. The frequency modulated M\"obius (FMM) model is defined as a parametric signal plus a gaussian noise, where the signal can be described as a single or a sum of waves. The FMM approach is flexible enough to describe a great variety of rhythmic patterns. The FMM package includes all required functions to fit and explore single and multi-wave FMM models, as well as a restricted version that allows equality constraints between parameters representing a priori knowledge about the shape to be included. Moreover, the FMM package can generate synthetic data and visualize the results of the fitting process. The potential of this methodology is illustrated with examples of such biological oscillations as the circadian rhythm in gene expression, the electrical activity of the heartbeat and neuronal activity.

2014], is widely used, in particular in chronobiology, with acceptable results when a sinusoidal shape response within a period is expected. Due to its widespread use, many software utilities are available. Particularly in R [R Core Team, 2020], the estimation of a COS model can be performed using cosinor [Sachs, 2014] and cosinor2 packages [Mutak, 2018]. In addition, other packages from widely differing areas of knowledge have specific functions for fitting COS models. Such is the case of, for example, the function CATCosinor in the CATkit package [Gierke et al., 2018], which implements tools for periodicity analysis; the function cosinor in the psych package [Revelle, 2018], dedicated to personality and psychological research; or the function cosinor contained in a very recent package, card [Shah, 2020], which is dedicated to the assessment of the regulation of cardiovascular physiology. Recently, it has also been implemented in other languages such as CosinorPy, a cosinor python package [Moskon, 2020]. The COS model is easy to use and interpret with symmetrical patterns. However, asymmetric shapes are not captured properly by COS. When the waveform is non-sinusoidal, the use of multiple components analysis to fit a model consisting of a sum of several periodical functions is recommended. However, the multicomponent FD models, developed to provide flexibility from COS, often require the use of a large number of components resulting in serious over-fitting issues.
In recent years, alternative methods, mostly nonparametric statistical methods, have been developed and used for analyzing rhythmicity, especially in biological data sets. Some very popular ones, such as the JTK_Cycle [Hughes et al., 2010], wrongly assume that any underlying rhythms have symmetric waveforms. Others, such as RAIN [Thaben and Westermark, 2014], designed to detect more diverse wave shapes including asymmetric patterns, are not focused on modeling but on detecting rhythmic behavior in sets of data. Thus, they are not useful to describe the underlying oscillatory phenomena. The proliferation of methodology in this field has been accompanied by software developments. This is the case, for example, of the DiscoRhythm R package [Carlucci et al., 2020], very recently available on Bioconductor (https://bioconductor.org/packages/DiscoRhythm) with a web interface based on the R Shiny platform [Chang et al., 2020]. This tool allows four popular approaches to be used, including the COS model and JTK_Cycle, to discover biological rhythmicity. Another recent example is the CircaCompare [Parsons et al., 2019], an R package implemented for modeling cosinusoidal curves by nonlinear regression and available from GitHub (https://github.com/RWParsons/circacompare). Also hosted on GitHub, we can find the LimoRhyde R package (https://github.com/hugheylab/limorhyde) for the differential analysis of rhythmic transcriptome data, based on fitting linear models [Singer and Hughey, 2019].
Motivated by the need for a flexible, interpretable and parametric methodology to fit rhythmic patterns, our research group recently proposed the frequency modulated Möbius (FMM) model [Rueda et al., 2019]. The FMM is an additive nonlinear parametric regression model capable of adapting to non-sinusoidal shapes and whose parameters are easily interpretable. The single component model has been shown to successfully fit data as diverse as circadian clock signals, hormonal levels data or light data from distant stars. In addition, for more complex oscillatory signals, a multicomponent model of order m, denoted as FMM m , which includes m single FMM components, can be used. This is, for example, the case for describing electrocardiography (ECG) signals. The FMM ecg signal, presented in Rueda et al. [2021b], is defined as the combination of five single FMM components. Another interesting area where the FMM approach has already shown its usefulness is in electrophysiological neuroscience. Specifically, we have proposed FMM methodology for modeling neuronal action potential (AP) curves, oscillating signals that measure the difference between the electrical potential inside and outside the cell [see Rueda et al., 2021c, Rodríguez-Collado and]. An FMM 2 model provides an accurate fitting for a single AP curve; whereas series of AP curves with similar repetitive spikes can be efficiently fitted by the FMM ST model, a restricted version of the multicomponent FMM model.
In this work we introduce the FMM package, programmed in R and available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=FMM. The package implements all required functions to fit and explore single and multicomponent FMM models, as well as a restricted multicomponent version. In addition, the FMM package provides functions to generate synthetic data and visualize the results of the fitted model. Furthermore, its use is illustrated in the aforementioned applications. The remainder of this paper is organized as follows: Section 2 provides a brief overview of both mono and multicomponent FMM models (2.1 and 2.2, respectively). The FMM m model with equality constraints is presented in Section 2.3. Section 3 is dedicated to the implementation details of the FMM package. Through simulated examples, Section 4 introduces the basic usage of the package, including the data generation and the fitting, as well as the visualization of the results. In Section 5, the FMM package performance is shown through three application areas governed by oscillatory systems: chronobiology, ECG and neuroscience. Finally, a summary is provided in Section 6.

Frequency modulated Möbius (FMM) model
All the methodological details that justify the mathematical formulation of the FMM models are given in Rueda et al. 2019. At the time point t, a single FMM wave is defined as W (t; υ) = A cos(φ(t; α, β, ω)) where υ = (A, α, β, ω) , A ∈ + represents the wave amplitude and, φ(t; α, β, ω) = β + 2 arctan(ω tan( t − α 2 )) (1) the wave phase. The phase angle φ of an FMM wave is defined using the Möbius link [see Downs andMardia, 2002, Kato et al., 2008] rather than the linear link function as in the COS model. The Möbius link provides much more flexibility to describe non-sinusoidal patterns. Without loss of generality, we assume that the time point t ∈ [0, 2π].
Each of the four parameters of an FMM wave characterizes some aspect of a rhythmic pattern. A describes the amplitude of the signal, while α, β and ω describe the wave phase. α ∈ [0, 2π] is a translation parameter and a wave location parameter in the real space, whereas β ∈ [0, 2π] and ω ∈ [0, 1] describe the wave shape. To be precise, assuming α = 0, the unimodal symmetric waves are characterized by values of β close to 0, π or 2π. When β = π 2 or β = 3π 2 , extreme asymmetric patterns are described. Moreover, a value of ω close to zero describes an extreme spiked wave and, as ω value increases, the pattern is increasingly smoother. When ω = 1, a sinusoidal wave is described and the FMM model matches the COS model where ϕ = β − α is the acrophase parameter.
Two important features of a wave are the peak and trough, defined as the highest and lowest points above and below the rest position, respectively. In many applications, the peak and trough times could be very useful tools to extract practical information of a wave, since they capture important aspects of the dynamics. These two interesting parameters can be directly derived from the main parameters of an FMM wave as, where t U and t L denote the peak and trough times, respectively.

Monocomponent FMM model
Let X(t i ), t 1 < t 2 < . . . < t n be the vector of observations. The monocomponent FMM model is defined as follows: where M ∈ is an intercept parameter describing the baseline level of the signal, W (t i ; υ) is an FMM wave, and it is assumed that the errors e(t i ) are independent and normally distributed with zero mean and a common variance σ 2 .

Estimation algorithm
A two-step algorithm to estimate monocomponent FMM model parameters is proposed. We now describe the substantial details of each stage of the algorithm.
Step 1: Initial parameter estimation. A two-way grid search over the choice of (α, ω) parameters is performed. For each pair of (α, ω) fixed values, the estimates for M , A and β are obtained by solving a least square problem as detailed below.
The model for a single FMM component can be written as: where t * i = α + 2 arctan(ω tan( ti−α 2 )), ϕ = β − α, and e(t i ) ∼ N (0, σ 2 ) for i = 1, . . . , n. Using trigonometric angle sum identity, the model can be rewritten as: where δ = A cos(ϕ), γ = −A sin(ϕ), z i = cos(t * i ) and w i = sin(t * i ). Since α and ω are fixed, the estimates for M , δ and γ are obtained by minimizing the residual sum of the squares (RSS), And the estimates for M , A and β are straightforward to derive as follows, The best combination of (α, ω) values, with the lowest RSS, is retained and the corresponding estimates are the initial parameter estimation values.
Step 2: Optimization. In the second step, the Nelder-Mead optimization method [Nelder and Mead, 1965] is used to obtain the final FMM parameter estimates. Nelder-Mead is a heuristic search method commonly used to minimize or maximize an objective function in a multidimensional space.

Multicomponent FMM model
A multicomponent FMM model of order m, denoted by FMM m , is defined as t 1 < t 2 < . . . < t n ; i = 1, . . . , n where W (t i ; υ J ), hereinafter denoted by W J (t i ), is the Jth FMM wave and,

Model adequacy
The goodness of fit of an FMM model is measured with the R 2 statistic that represents the proportion of the variance explained by the model out of the total variance, that is: whereX(t i ) represents the fitted value at t i , i = 1, . . . , n.

Estimation algorithm
An iterative backfitting algorithm is proposed to derive estimates for the FMM parameters. LetŴ and fit a monocomponent FMM model to r Repeat the backfitting step until the stopping criterion is reached. The stopping criterion is defined as the difference between the explained variability in two consecutive iterations: Equation 11) is the proportion of variance explained by the model in the kth iteration and C a constant.

Restricted multicomponent FMM model
Modeling signals with repetitive shape-similar waves can be very useful in some applications [see Rodríguez-Collado and Rueda, 2021b]. In order to obtain more efficient estimators, equality constraints are imposed on the β and ω parameters of an FMM m model. In particular, we add d blocks of restrictions on β parameters: . . . β m d−1 +1 = . . . = β m d and/or d blocks of restrictions on ω parameters: The parameter estimation problem is solved by an adaptation of the standard procedure.

FMM m estimation algorithm with restrictions on the β parameters
Given the unrestricted estimates obtained in step 3, the estimates for β 1 , β m1+1 , . . . , β m d−1 +1 under equality restrictions (Equation 14) are computed as follows: Then, the algorithm continues to the next step.

FMM m estimation algorithm with restrictions on the ω parameters
A block is defined as the set of FMM waves with common ω parameter. When constraints for the ω parameters are incorporated, the grid search for the different ω values is outside the backfitting loops. When the number of blocks is large, the estimation procedure can be computationally unaffordable. In order to reduce the execution time, a two-nested backfitting algorithm is proposed. In the outer backfitting loop, a block is fitted. In the inner loop, the FMM waves belonging to the same block are estimated. This procedure generates a close to optimal solution and is a less computationally expensive alternative.

FMM package: Implementation details
The FMM code makes use of the doParallel package [Corporation and Weston, 2020] to embed parallelization for the restricted version of the fitting process. Several utilities from the ggplot2 [Wickham, 2016] and RColorBrewer [Neuwirth, 2014] packages are occasionally necessary for the visualization of the fitted models.
The implementation of FMM is divided into four main functionalities described from Sections 3.1 to 3.4: the fitting of the FMM models, the new S4 object of class 'FMM', the graphical visualization of the fittings and the simulation of synthetic data.
Some general details about the functions contained in the FMM package are shown in Table 1.

Fitting an FMM model
An FMM model can be fitted using the main function fitFMM(). The description and default values of its inputs arguments are shown in Table 2.

extractWaves(objFMM)
Extracts individual contribution to the fitted values of each FMM wave. Standard methods for objects of class 'FMM' summary(), show(), coef(), fitted() Table 1: Summary of the fitting, utility functions and standard methods implemented in FMM package.
The fitting function fitFMM() requires the vData input argument, which contains the data to be fitted. Two other arguments can be used to control a basic fitting: timePoints, which contains the specific time points of the single period; and nback, with the number of FMM components to be fitted. For some applications, such as the study of circadian rhythms, data are collected over multiple periods. This information is received by the fitFMM() function through the input argument nPeriods. When nPeriods > 1, the FMM fitting is carried out by averaging the data collected at each time point across all considered periods. There are three key issues in the fitting process: the grid search of the pair (α, ω) to solve the estimation problem of a single FMM wave, the backfitting algorithm used for the estimation of the multicomponent models, and the incorporation of restrictions on β and ω parameters. Each of these issues is controlled by several arguments described below.

Argument
• Grid search of the pair (α, ω). The lengthAlphaGrid and lengthOmegaGrid arguments are used to set the grid resolution by specifying the number of equally spaced α and ω values. Thus, the objetive function will be evaluated a total number of (lengthAlphaGrid) × (lengthOmegaGrid) times, so when both arguments are large, the computational demand can be high. By reducing the size of the sequences of the α and ω parameters, the algorithm will be computationally more efficient. However, it may fail to obtain an accurate estimation if the grid resolution is too sparse. An implemented option to fine-tune the estimation of the parameters is to repeat the fitting process a numReps of times, in such a way that, at each repetition, a new two-dimensional grid of (α, ω) points is created around the previous estimates. • Backfitting algorithm. The argument maxiter sets the maximum number of backfitting iterations. Through the argument stopFunction, it is possible to set a stopping criterion. Two criteria have been implemented as stop functions in this package. When stopFunction = alwaysFalse, maxiter iterations will be forced. If stopFunction = R2(), the algorithm will be stopped when the difference between the explained variability in two consecutive iterations is less than a value pre-specified in the difMax argument.
• Restrictions. The arguments betaRestrictions and omegaRestrictions set the equality constraints for the β and ω parameters, respectively. For the unrestricted case, betaRestrictions = 1 : nback and omegaRestrictions = 1 : nback. To add restrictions, 'integer' vectors of length m can be passed to these arguments, so that positions with the same numeric value correspond to FMM waves whose parameters, β and/or ω, are forced to be equal. Since restricted fitting can be computationally intensive, especially for constraints on ω parameters, a parallel processing implementation can be used when the argument parallelize = TRUE.

Object of class 'FMM'
The fitFMM() function outputs an S4 object of class 'FMM' which contains the slots presented in Table 3. A m-element 'numeric' vector specifying the explained variance by each of the fitted FMM components. nIter A 'numeric' value containing the number of iterations of the backfitting algorithm. Table 3: Summary of the slots of the S4 object of class 'FMM' resulting from fitting an FMM model with m components.
The standard methods implemented for the class 'FMM' include the functions summary(), show(), coef() and fitted(). These methods display relevant information of the FMM fitting, and provide the estimated parameters and fitted values. In addition, two more specific functions have been implemented. Through the extractWaves() function, the individual contribution of each FMM wave to the fitted values can be extracted. Finally, the location of the peak and trough of each FMM wave, as well as the value of the signal at these time points, can be estimated using the getFMMPeaks() function. The required argument of all these methods and functions is an object of the class 'FMM'. Particularly, getFMMPeaks() has an optional argument: timePointsIn2pi, that forces the peak and trough locations to be returned into the interval from 0 to 2π when it is TRUE.

Plotting FMM models
The FMM package includes the function plotFMM() to visualize the results of an FMM fit. The arguments of this function are summarized in Table 4.
An object of class 'FMM' can be plotted in two ways (see Figure 1). The default graphical representation will be a plot on which original data (as points) and the fitted model (as a line) are plotted together (left panel in Figure 1). The other possible representation is a component plot for displaying each centered FMM wave separately (right panel in Figure 1). Set the boolean argument components = TRUE to show a component plot. When legendInComponentsPlot = TRUE, a legend appears at the bottom of the component plot to indicate the represented waves. The argument textExtra allows an extra text to be added to the title of both graphical representations.
As mentioned above, in some cases, data are collected from different periods. All periods can be displayed simultaneously on the default plot using plotAlongPeriods = TRUE. For the component plot, this argument is ignored. The argument use_ggplot2 provides a choice between building the plot using base R graphics or ggplot2 packages. By default, the graphics package is used. When use_ggplot2 = TRUE, a more aesthetic and customizable plot is created using the ggplot2 package.

Simulating data from an FMM model
Data from an FMM model can be easily simulated using the function generateFMM() of the package FMM. All input arguments of this function are shown in Table 5 The main arguments of this function are M, A, alpha, beta and omega, whereby the values of the FMM model parameters are passed to the function. All these arguments are 'numeric' vectors of length m, except M, which has length 1. Longer and smaller vectors will be truncated or replicated as appropriate.
By default, the data will be simulated at a sequence of 100 equally spaced time points from 0 to 2π. The arguments from, to and length.out control such sequences. The sequence can also be manually set using the argument timePoints, in which case from, to and length.out will be ignored.
The user can add a gaussian noise by argument sigmaNoise. A positive 'numeric' value sets the corresponding standard deviation of the gaussian noise to be added. To create the normally distributed noise, the rnorm() function is used.
The arguments plot and outvalues, both boolean values, determine the output of the generateFMM() function. When outvalues = TRUE, a 'list' with input parameters, time points and simulated data is returned. These elements are named input, t and y, respectively. In addition, a scatter plot of y against t can be drawn by setting plot = TRUE.

Basic usage of the FMM package
The two examples bellow, based on FMM synthetic data, illustrate the basic uses and capabilities of the functions implemented in the FMM package. The first example is based on a simulated data set from a multicomponent FMM model and illustrates the fitting process and the presentation of its results. The second is devoted to the restricted version of FMM m models and serves to exemplify the control of the grid search related arguments to reduce the execution time.

Example 1: A multicomponent FMM data set
The first simulated data set consists of 100 observations from an FMM 2 model with intercept parameter M = 0 and two FMM waves defined by the parameters υ 1 = (A 1 = 2, α 1 = 1.5, β 1 = 0.2, ω 1 = 0.1) and υ 2 = (A 2 = 2, α 2 = 3.4, β 2 = 2.3, ω 2 = 0.2) , respectively. The standard deviation of the error term is set at σ = 0.3. We use the function generateFMM() to simulate this data set. A set.seed() statement is used to guarantee the reproducibility of the results. When the argument showProgress = TRUE (its default value), a progress indicator and information about the stopping criterion of the backfitting algorithm will be displayed on the console. By default, maxiter iterations will be forced, but we can use the argument stopFunction to modify the criteria. For example, to continue the search until there is an improvement, in terms of explained variability, of less than 1%, the R2() function can be used specifying the difMax = 0.01 argument. The FMM wave parameter estimates, as well as the peak and trough times, together with the signal values at those times, are presented in tabular form, where each row corresponds to a component and each column to an FMM wave parameter. As part of the summary, a brief description of the residuals, the proportion of variance explained by each FMM component and by the global model are also shown. From the above summary results, we can see that the variance explained by the fitted model is 96.85% and that the FMM waves are labelled in decreasing order according to the R 2 value: the explained variance is 69.07% and 27.78% by "FMM wave 1" and "FMM wave 2", respectively. The summary output can be assigned to an object to get a 'list' of all the displayed results.
Other options to return the results are the functions coef(), getFMMPeaks() and resid(). The first two return a 'list' similar to those obtained with summary(). The resid() method can be used to obtain the complete residuals vector. In addition, the fitted values can be extracted by the function fitted(), which returns a 'data.frame' with two columns: time points and fitted values.

R> f i t t e d . fmm2 <− f i t t e d ( f i t . fmm2 ) R> h e a d ( f i t t e d . fmm2 )
t i m e P o i n t s f i t t e d 1 0 . 0 0 0 0 0 0 0 0 1 . 2 5 1 9 3 7 2 0 . 0 6 3 4 6 6 5 2 1 . 7 2 6 5 7 0 3 0 . 1 2 6 9 3 3 0 4 2 . 2 6 9 9 2 8 4 0 . 1 9 0 3 9 9 5 5 2 . 8 2 6 4 4 4 5 0 . 2 5 3 8 6 6 0 7 3 . 3 1 9 6 8 3 6 0 . 3 1 7 3 3 2 5 9 3 . 6 8 1 0 2 4 The code given below creates Figure 1. R> l i b r a r y ( " R C o l o r B r e w e r " ) R> l i b r a r y ( " g g p l o t 2 " ) R> l i b r a r y ( " g r i d E x t r a " ) R> t i t l e T e x t <− " Example 1 : two s i m u l a t e d FMM waves " R> defaultFMM2 <− plotFMM ( f i t . fmm2 , u s e _ g g p l o t 2 = TRUE , t e x t E x t r a = t i t l e T e x t ) + theme ( p l o t . m a r g i n = u n i t ( c ( 1 , 0 . 2 5 , 1 . 3 , 1 ) , "cm " ) ) + y l i m ( − 5 , 5 ) R> compFMM2 <− plotFMM ( f i t . fmm2 , c o m p o n e n t s = TRUE , u s e _ g g p l o t 2 = TRUE , t e x t E x t r a = t i t l e T e x t ) + theme ( p l o t . m a r g i n = u n i t ( c ( 1 , 0 . 2 5 , 0 , 1 ) , "cm " ) ) + y l i m ( − 5 , 5 ) + s c a l e _ c o l o r _ m a n u a l ( v a l u e s = b r e w e r . p a l ( " S e t 1 " , n = 8 ) [ 3 : 4 ] ) R> g r i d . a r r a n g e ( defaultFMM2 , compFMM2 , nrow = 1 ) The FMM plots can be written in the R graphics or ggplot2 packages. In this example, we use use_ggplot2 = TRUE to build Figure 1 based on ggplot2. The use of ggplot2 makes it easier to customize our plots and modify features, such as scales, margins, axes, etc. In Figure 1, the two possible FMM plots are arranged via the grid.arrange function of the gridExtra package [Auguie, 2017].

Example 1: Chronobiology
Chronobiology studies ubiquitous daily variations found in nature and in many aspects of the physiology of human beings, such as blood pressure or hormone levels [Mermet et al., 2017]. These phenomena commonly display signals with oscillatory patterns that repeat every 24 hours, usually known as circadian rhythms. In particular, circadian gene expression data have been deeply analyzed in the literature as they regulate the vast majority of molecular rhythms involved in diverse biochemical and cellular functions, see among others Zhang et al. [2014], Cornelissen [2014] and Larriba et al. [2020].
The FMM package includes a data set called mouseGeneExp that contains expression data of the Iqgap2 gene from mouse liver. The liver circadian database is widely extended in chronobiology since the liver is a highly rhythmic organ with moderate levels of noise [Anafi et al., 2017, Larriba et al., 2018. The complete database is freely available at NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/), with GEO accession number GSE11923. Gene expression values are given along 48 hours with a sampling frequency of 1hour/2days. Hence, data are collected along two periods, and an FMM 1 model is fitted to the Iqgap2 average expressed values as follows: R> d a t a ( " mouseGeneExp " , p a c k a g e = "FMM" ) R> f i t G e n e <− fitFMM (  The behavior of the FMM versus COS model to describe this asymmetric pattern has been compared in terms of R 2 . The FMM model clearly outperforms the COS one with an R 2 of 0.8752 and 0.2835, respectively. In addition, a difference of 4.73 hours in peak time estimation between both models is observed, the FMM peak estimate being much more reliable, as is shown in Figure 3.

Example 2: Electrocardiography
ECG records the periodic electrical activity of the heart. This activity represents the contraction and relaxation of the atria and ventricle, processes related to the crests and troughs of the ECG waveform. Heartbeats are decomposed into five fundamental waves, labelled as P , Q, R, S and T , corresponding to the different phases of the heart's electric activity. The main features used in medical practice for cardiovascular pathology diagnosis are related to the location and amplitudes of these waves, and, of them, those labeled as P , R and T are of particular interest [Bayes de Luna, 2007]. Standard ECG signals are registered using twelve leads, calculated from different electrode locations, Lead II being the reference one.
The FMM package includes the analysis of a typical ECG beat from the QT database [Laguna et al., 1997]. This recording, from the subject sel100, belongs to the Normal category, regarding Physionet's pathology classification [Goldberger et al., 2000]. The data illustrate the voltage of the heart's electric activity, measured in mV , along the beat with a sampling frequency of 250Hz. Specifically, the ECG signal from Lead II in the fifth of the thirty annotated beats is analysed. Recordings are publicly available on (http://www.physionet.org). Data are saved as ecgData in the package. For an ECG beat, an FMM ecg , a fifth order multicomponent FMM model can be fitted with the instruction: R> d a t a ( " e c g D a t a " , p a c k a g e = "FMM" ) R> f i t E c g <− fitFMM (  It is worth noting that the FMM package not only provides ECG signal-fitting (the left hand panel in Figure 4), but it also does wave decomposition and fiducial mark annotations on the desired waves (the right hand panel in Figure 4). It is clearly visible how the specific shapes of the five main waves contribute to drawing and explaining the Lead II ECG waveform from the Normal morphology. See Rueda et al. [2021b] for a complete review of FMM ecg .

Single AP curve
The study of the electrophysiological activity of neurons is one of the main research branches in neuroscience. The AP curves are oscillatory signals that serve as basic information units between neurons. They measure the electrical potential difference between inside and outside the cell due to an external stimulus. Gerstner et al. [2014] can serve as a basic reference for electrophysiological neuroscience. Recently, the shape and other features of the AP have been used in problems such as spike sorting [Rácz et al., 2020, Souza et al., 2019, Caro-Martín et al., 2018 or neuronal cell type classification [Teeter et al., 2018, Gouwens et al., 2019, Mosher et al., 2020, Rodríguez-Collado and Rueda, 2021a.
The package includes an example of a neuronal AP. The data were simulated with the renowned Hodgkin-Huxley model, first presented in Hodgkin and Huxley [1952], which is defined as a system of ordinary differential equations and has been used in a wide array of applications, as it successfully describes the neuronal activity in various organisms. The goodness of fit of the FMM 2 model can be ascertained in Figure 5. For comparison purposes, an FD model has been fitted with the same number of degrees of freedom. While the FD attains an R 2 = 0.3926, the FMM model achieves a better fit with R 2 = 0.9669.  Figure 5: Neuronal AP simulated with the Hodgkin-Huxley model (parameters: C = 1, g N a = 260, g K = 30, g L = 0.31, V K = −12, V N a = 115, V L = 10.6, a n = 1.15, b n = 0.85, a m = 0.9, b m = 1.3, a h = 1, b h = 1 and applying a current of 12µA for 1 millisecond) and the estimated FMM 2 signal in red. An FD model of the same number of degree of freedom has been fitted and plotted in blue.

AP train
Multiple AP curves, denominated spike or AP train, are usually observed as the response to a stimulus. Various models, such as the widely used leaky-and-fire models [Lynch and Houghton, 2015], cut the signal into segments, each one containing an AP curve. Some authors suggest cutting the signal into even segments [Gerstner et al., 2014]. However, the length of the segments turns out to be significantly different between different types of neurons, as explained in Teeter et al. [2018], and unequal data segments can lessen the utility of some approaches. An important aspect to take into account is that the shape of the APs in the spike train is considered to be similar and, consequently, a restricted FMM model can accurately fit the entire signal.
The FMM package includes the data of a spike train composed of three AP curves. The proposed model for use with these data is an FMM ST model, as defined in Rodríguez-Collado and . Each AP is modeled by two components. The β and ω parameters are constrained between AP curves. The code below fits the model.
The package allows both single and multicomponent FMM models to be estimated. In order to provide greater flexibility, equality constraints for shape parameters have also been implemented. In addition, graphical representations of the fitted models and the possibility of generating synthetic data are available. The functionality of the package has been illustrated by simulated data and also by real examples from different areas of application related to present-day biological problems.
Possible future extensions of the FMM package include the implementation of additional restrictions to suit the model to other real signals; the possibility to include weights that determine how much each observation influences the parameter estimates; and the improvement of the computational efficiency. Currently, one of the most important limitations of the stable version of the package (version 0.1.2) is that the fitting process, especially of the restricted model, can be a time-consuming process. Nevertheless, the execution times of the development version available on the GitHub repository (https://github.com/alexARC26/FMM) have already been significantly reduced, in many cases by an order of 10 times or more.