The highest values in terms of sensitivity were obtained by NB in WC, BMI, and VAI (66.9%, 69.7%, and 70.5%, respectively), and by RF in WHtR, BRI, and CUN-BAE (68.7%, 69.4%, and 68.3%, respectively). Regarding specificity in all the anthropometric parameters, the highest values were obtained by LR. Sensitivity for all three algorithms and for all anthropometric parameters showed similar results, with values ranging from 0.650 to 0.705. The highest values in terms of sensitivity for the anthropometric parameters VAI, BMI, and WC (70.5%, 69.7%, and 66.9%, respectively) were obtained with NB. The highest values in terms of sensitivity for BRI, WHtR, and CUN-BAE (69.4%, 68.7%, and 68.3%, respectively) were achieved with RF. Regarding specificity in all the anthropometric parameters, the highest values were obtained by LR, ranging from 0.488 to 0.528. RF and NB presented very low values in specificity in all anthropometric parameters, with the exception of the NB algorithm for the anthropometric parameter VAI.

#### **4. Discussion and Conclusions**

Through a logistic regression model, we explored the capacity of anthropometric parameters to predict vitamin D deficiency, and we concluded that the behavior of these parameters differed according to sex; thus, in males, WC, WHtR, VAI, and BRI were associated with low levels of vitamin D, whereas in females, CUN-BAE was associated with low vitamin D.

To the best of our knowledge, this was the first study to use ML algorithms for the detection of vitamin D and investigate its predictability using anthropometric parameters.

Vitamin D deficiency (<20 ng/mL) in our study affected 34.7% of the participants. These results were in line with the systematic review conducted by Manios et al. [56], who reported that in Southern European countries, over one-third of the population had vitamin D levels <20 ng/mL, and 10% of the population had values <10 ng/mL. A cross-sectional, retrospective study with 21,490 patients (74.3% females) aged between 14 and 105 years who had used primary healthcare in La Rioja (Spain) showed that the mean levels of 25(OH)D were 18.3 (SD, 11.6) ng/mL in the entire sample [57], with males presenting lower values than females (17.6 vs. 18.5 ng/mL, *p* < 0.001). However, our study did not obtain significant differences according to sex, and the mean levels were above 25.56 (SD, 19.30), likely due to the differences in age, sex, and associated diseases between the study populations.

The relationship between low concentrations of 25(OH)D and obesity has been previously reported [58,59]. Jääskeläine et al. [60] used data from the 2000–2011 Health Survey and suggested that vitamin D deficiency may be a risk factor for abdominal obesity among males, but not among females. These results were in line with those of our study; in males, in the logistic regression analysis, WC and BMI were associated with vitamin D deficiency. Nevertheless, these results disagreed with those of Cătoi et al. [61], who explored the complex relationship between the levels of 25(OH)D and overweight/obesity, insulin resistance, systemic inflammation, and oxidative stress, revealing that overweight and an increasing degree of obesity were not significantly associated with a decrease in the levels of 25(OH)D.

A Danish study with 4909 children and adolescents in the Danish Childhood Obesity Biobank (2860 females) found that vitamin D deficiency was common among Danish children and adolescents with obesity. Our results were in line with those of that study, as the individuals with vitamin D deficiency showed greater BMI, WHtR, and frequency in males, which indicated that the degree of obesity was independently associated with lower serum concentrations of 25(OH)D [62]. Moreover, it is known that obese people need higher vitamin D loading doses to reach the same amount of serum 25-hydroxyvitamin D as people with normal body weight [5]. However, not all studies have found this association. The results of Pereira Santos et al. [4] indicated that overweight and obese individuals in different age groups have a similar probability of presenting with vitamin D deficiency. These discrepancies could be due to the origins of the different study populations.

In a study conducted with young Italian females, Adami et al. [15] found that the main determinants of vitamin D deficiency were an increase in BMI and exposure to sunlight. In a more recent study performed to identify the best combination of predictors for the serum concentration of vitamin D in adults aged between 18 and 70 years old, the multivariate linear regression model included age, sex, BMI, sunlight exposure in the previous week and during the month of blood sample collection, skin phototype, job position, smoking status, physical activity, latitude, and administration of vitamin D supplements in the previous year [63]. However, these results differed from some of our results, in that BMI was not found to be significant in females; we considered that these differences could be due to differences in age ranges and ethnicities.

The relationship we found between vitamin D deficiency and VAI was in disagreement with the study by Izadi et al. [64], who analyzed a sample of 57 males and 26 females with nonalcoholic fatty liver disease (NAFLD) and despite controlling for age and sex, they found a reverse association between VAI and vitamin D levels. Nevertheless, a study by Zubiaga et al. [65] reported that body fat percentage (BFP), when calculated with CUN-BAE as a predictive marker of cardiovascular risk in patients with morbid obesity before and after being subjected to vertical gastrectomy (VG), was significantly correlated with three biochemical factors associated with greater cardiovascular risk (i.e., cortisol, vitamin D, and TG/HDL-C ratio), which was in line with our results, as we only found an association between CUN-BAE and vitamin D in females.

Several studies have analyzed the association (i.e., predictive ability) of 25-hydroxyvitamin D concentrations with various health problems, using different ML techniques. Luo et al. [66] analyzed whether 25-hydroxyvitamin D deficiency was associated with an increased incidence of COVID-19 and disease severity using multivariable logistic regression techniques in order to propose a predictive model. These authors observed that subjects with COVID-19 had lower 25-hydroxyvitamin D concentrations as compared to the controls without COVID-19, and that 25-hydroxyvitamin D deficiency influenced both hospitalization rates and the severity of COVID-19 in the Chinese subjects. Deschasaux et al. [67] used logistic regression in 1557 middle-aged adults without prior 25-hydroxyvitamin D treatment to develop a scale to predict 25-hydroxyvitamin D deficiency and identify adults at risk of deficiency. This scale indicated that in subjects with scores of ≥7 points, 70% were deficient in 25-hydroxyvitamin D, and when the score was >9, 80% were deficient in 25-hydroxyvitamin D, with a sensitivity of 0.67 and a specificity of 0.63. Therefore, the application of this scale could avoid unjustified 25-hydroxyvitamin D supplementation and unnecessary blood tests. Garcia-Carretero et al. [37] analyzed 1002 hypertensive patients to establish predictive models to identify patients unlikely to have 25-hydroxyvitamin D deficiency or to undergo plasma 25-hydroxyvitamin D concentration measurements. To do so, they used the classifiers logistic regression, support vector machine (SVM), RF, NB, and extreme gradient boosting to calculate classification accuracy, sensitivity, specificity, and predictive values to assess the performance of each method. These authors found that the radial kernel, SVM-based method performed better than the other algorithms in terms of sensitivity (98%), negative predictive value (71%), and classification accuracy (73%). Therefore, they concluded that the combination of a feature-selection method such as elastic regularization, as well as a classification approach, produced well-fitted models. This combined approach allowed them to develop a prediction model with high sensitivity but low specificity, which was consistent with the results of our study, to identify the population that could benefit from a laboratory determination of serum 25-hydroxyvitamin D levels. Guo et al. [38] analyzed MLR and RBF-SVR techniques in 594 Caucasian adults to develop a score for predicting serum 25-hydroxyvitamin D concentration. The best results were found using the RBF-SVR model, which provided a better prediction of serum 25-hydroxyvitamin D concentrations and vitamin D deficiency, as compared to an MLR model. Lopes et al. [68] analyzed 908 community-dwelling older people using logistic regression to propose a model for detecting 25-hydroxyvitamin D deficiency. The model was able to identify older people at risk of 25-hydroxyvitamin D deficiency with a sensitivity of 55.9%, a specificity of 72.3%, and an ROC area of 0.685. These authors suggested that a clinical use of these parameters could help to identify and design appropriate public health interventions. Finally, Sohl et al. [69], in a longitudinal aging study in Amsterdam with 1509 subjects, developed a risk profile based on backward logistic regression to identify older people at high risk of 25-hydroxyvitamin D deficiency. In this study, two total risk scores were developed that included either 10 or 13 variables that were capable of predicting serum 25(OH)D concentrations of less than 0.50 and 0.30 nmol/L, respectively. This scale may be useful in clinical practice to identify individuals at risk of 25-hydroxyvitamin D deficiency.

Given the association between 25-hydroxyvitamin D deficiency and obesity, and that healthcare professionals deal with many variables that can influence health problems, ML algorithms have potential clinical application.

Recent research has further explored linking these techniques to provide hybrid ML algorithms. Therefore, the approach of this paper, which indicated that the predictive ability of different anthropometric parameters differed according to sex, could be useful in future research. However, additional studies are needed to confirm our results.

#### *4.1. Limitations of the Study*

This study had several limitations. Firstly, this was a cross-sectional study, which hindered the establishment of causal relationships between vitamin D levels and anthropometric parameters. Secondly, there may have been confounding variables that were not considered in this study. Lastly, the number of patients with vitamin D deficiency in our study was unbalanced as compared to the number of people with normal levels.

#### *4.2. Conclusions*

The capacity of anthropometric parameters to predict vitamin D deficiency differed according to sex; thus, WC, BMI, WHtR, VAI, and BRI were useful predictors in males, while CUN-BAE was more useful in females. In all the anthropometric parameters, the LR model presented the highest values in terms of specificity to predict vitamin D deficiency. The NB approach of ML showed the best area under the curve in WC, BMI, WHtR, and BRI, whereas the LR model did so for VAI and CUN-BAE.

**Author Contributions:** Conceptualization, C.P.-A., M.A.G.-M., L.G.-O. and E.R.-S.; methodology, M.G.-S., L.G.-S., M.A.G.-M., E.R.-S. and B.S.S.; data analysis, C.P.-A., M.A.G.-M. and L.G.-O.; writing review and editing, all authors; supervision, C.P.-A.; funding acquisition, M.A.G.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** The project was funded by the Institute of Health Carlos III of the Spanish Ministry of Science, Innovation, and Universities through Red de investigación en cronicidad, atención primaria y promoción de la salud (RD21/0016/0010), and cofinanced by the European Union Health Institute/European Regional Development Fund (ERDF), the Autonomous Government of Castilla and León (GRS 1193/B/15, GRS 1821/B/18), and intensification of a research program (INT/M/08/19, INT/M/9/19/INT/M/14/19).

**Institutional Review Board Statement:** The study was approved by the Drug Research Ethics Committee of the Health Area of Salamanca on 4 May 2015 under registry number PI15/01039. All procedures conducted with human participants complied with the ethical rules of the institutional and/or national research committee and with the 2013 Declaration of Helsinki [70]. All participants signed an informed consent form before participating in this study.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare that they have no known competing financial interest or personal relationships that could have influenced the work reported in this paper.

#### **References**

### *Article* **Effect of Probability Distribution of the Response Variable in Optimal Experimental Design with Applications in Medicine †**

**Sergio Pozuelo-Campos \*,‡, Víctor Casero-Alonso ‡ and Mariano Amo-Salas ‡**

Department of Mathematics, University of Castilla-La Mancha, 13071 Ciudad Real, Spain<EMAIL_ADDRESS>(V.C.-A<EMAIL_ADDRESS>(M.A.-S.)

**\*** Correspondence<EMAIL_ADDRESS>
† This paper is an extended version of a published conference paper as a part of the proceedings of the 35th International Workshop on Statistical Modeling (IWSM), Bilbao, Spain, 19–24 July 2020.

‡ These authors contributed equally to this work.

**Abstract:** In optimal experimental design theory it is usually assumed that the response variable follows a normal distribution with constant variance. However, some works assume other probability distributions based on additional information or practitioner's prior experience. The main goal of this paper is to study the effect, in terms of efficiency, when misspecification in the probability distribution of the response variable occurs. The elemental information matrix, which includes information on the probability distribution of the response variable, provides a generalized Fisher information matrix. This study is performed from a practical perspective, comparing a normal distribution with the Poisson or gamma distribution. First, analytical results are obtained, including results for the linear quadratic model, and these are applied to some real illustrative examples. The nonlinear 4-parameter Hill model is next considered to study the influence of misspecification in a dose-response model. This analysis shows the behavior of the efficiency of the designs obtained in the presence of misspecification, by assuming heteroscedastic normal distributions with respect to the D-optimal designs for the gamma, or Poisson, distribution, as the true one.

**Keywords:** elemental information matrix; gamma distribution; poisson distribution; D-optimization; misspecification

#### **1. Introduction**

To obtain optimal designs, it is common to assume a homoscedastic normal distribution of the response variable and under this assumption there is vast literature focused mainly on nonlinear models. However, there are also papers that use probability distributions different from a normal distribution [1–7]. At this point, it is important to remember that the probability distribution of the response variable is assumed, on many occasions, from the nature of the experiment to be performed. However, there are usually no prior observations to allow this assumption to be checked.

There are very few available references that set out a general framework for optimal experimental design for any probability distribution of the response variable. Ref. [8] present a method to compute the D-optimal designs for Generalized Linear Models with a binary response allowing uncertainty in the link function, ref. [9] study the Generalized Linear Model from the perspective of optimal experimental design, ref. [10] present the "elemental information matrix" for different probability distributions, and [11] compute optimal designs based on the maximum quasi-likelihood estimator to avoid the misspecification in the probability distribution of the response. The aim of this paper is to analyze the effect of misspecification in the probability distribution in optimal design. In other words, it allows those cases to be identified in which it is important to pay special attention to the assumed probability distribution. In this study, apart from theoretical results, real

**Citation:** Pozuelo-Campos, S.; Casero-Alonso, V.; Amo-Salas, M. Effect of Probability Distribution of the Response Variable in Optimal Experimental Design with Applications in Medicine. *Mathematics* **2021**, *9*, 1010. https:// doi.org/10.3390/math9091010

Academic Editor: Lev Klebanov

Received: 18 March 2021 Accepted: 27 April 2021 Published: 29 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

applications involving the linear quadratic model and a dose-response model are considered. For the latter, we focus on the well-known Hill model, widely used to describe dependence between the concentration of a substance and a variety of responses in biochemistry, physiology or pharmacology. From the point of view of optimal experimental design, this model is studied in many papers [12–15]. Specifically, ref. [13] study the effect of some drugs which inhibit the growth of tumor cells providing D-optimal designs under the assumption of the response variable follows a heteroscedastic normal distribution with a given structure for the variance.

The article is organized as follows. Section 2 introduces the model used and the theory of optimal experimental design. Section 3 presents the structure of the variance of the heteroscedastic normal distribution and proves a general theoretical result. Section 4 focuses on the linear quadratic model and provides some theoretical results for gamma or Poisson distributions. This section also shows applications of these results to real examples found in the literature. Finally, the 4-parameter Hill model is studied in Section 5. Assuming the heteroscedastic normal distribution, as in [13], an efficiency analysis is performed, considering the Poisson, or gamma, as the true probability distribution. A sensitivity analysis with respect to a parameter of the variance structure is also performed. The paper concludes with a summary and conclusions section.

#### **2. Model and Optimal Experimental Design**

The model of interest to the practitioner is expressed in a general way as

$$E[y] = \mathfrak{g}^{-1}(\eta(\mathfrak{x}; \mathfrak{e})),\tag{1}$$

where *y* is the response variable, following a probability distribution with pdf *d*(*y*; *ρ*), where *ρ* is the vector of parameters of the assumed distribution, *η*(*x*; *θ*) is the regression function (linear or nonlinear in the parameters), *x* is the vector of controllable variables and *θ* the vector of unknown parameters that must be estimated. Lastly, *g* is the link function relating the regression function to the mathematical expectation of the response. Ref. [16] carry out an in-depth study of the link function and Generalized Linear Models. In line with these authors, this paper considers the canonical link function for the probability distributions involved in the study, as it guarantees that the maximum likelihood estimators of the model parameters, *θ*ˆ, are sufficient.

An exact design of size *n* is defined as a set of values of the explanatory variables, *x*1, ... , *xn*, in which some may be repeated. These values belong to a compact set called design space <sup>X</sup> , which is usually a subset of <sup>R</sup>*N*. However, the real applications, examples and results in this study consider the one-dimensional case. Assuming that only *q* of these values are distinct, we may consider the set *x*1, ... , *xq* and associate with it a probability measure defined by *w*1, ... , *wq*, where each *wi* represents the proportion of experiments carried out under the condition *xi*. This suggests a more general definition of approximate design as a probability measure *ξ* over the design space X :

$$\xi = \begin{Bmatrix} \mathfrak{x}\_1 & \dots & \mathfrak{x}\_q \\ w\_1 & \dots & w\_q \end{Bmatrix} \in \Xi, \quad \sum\_{i=1}^q w\_i = 1,$$

where *<sup>ξ</sup>*(*xi*) = *wi* and <sup>Ξ</sup> represents the set of all approximate designs.

The scenario studied in this work is the estimation of a single parameter of the probability distribution of the response, with the rest being fixed. Thus, the elemental information matrix (EIM), introduced by [10], is scalar and is defined as

$$\nu(\eta(\mathbf{x};\boldsymbol{\theta})) = -E\left[\frac{\partial^2 \log d(\underline{y}; \eta(\mathbf{x}; \boldsymbol{\theta}))}{\partial \eta(\mathbf{x}; \boldsymbol{\theta})^2}\right],\tag{2}$$

which contains information about the probability distribution of the response variable *y*, given by the pdf *d*(*y*; *ρ*). The relationship between the parameters to estimate, *ρ*, of the probability distribution and the regression function *η*(*x*; *θ*) is established by the link function, *g*, shown in (1). Table 1 sets out the canonical link function, the mathematical expectation of the response variable as a function of *η*(*x*; *θ*) and the EIM for the probability distributions used in this paper, some of which are derived in Section 3.

**Table 1.** Density function, link function, expectation of the response variable as a function of *η* = *η*(*x*; *θ*) and the EIM for the probability distributions used in this paper.

The single-point information matrix in *x* ∈ X is given by

$$I(\mathbf{x};\boldsymbol{\theta}) = -E\left[\frac{\partial^2 \log d(\mathbf{y}; \boldsymbol{\eta}(\mathbf{x}; \boldsymbol{\theta}))}{\partial \theta\_i \partial \theta\_j}\right] = \nu(\boldsymbol{\eta}(\mathbf{x}; \boldsymbol{\theta})) f^T(\mathbf{x}; \boldsymbol{\theta}) f(\mathbf{x}; \boldsymbol{\theta}), \ \forall i, j = 1, \dots, m\_i$$

where *ν*(*η*(*x*; *θ*)) is the EIM defined in (2) and

$$f^T(\mathbf{x}; \boldsymbol{\theta}) = \frac{\partial \eta(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}}.$$

Finally, the Fisher information matrix (FIM) is defined for the approximate design with probability measure *ξ* as

$$\mathbf{M}(\xi;\theta) = \int\_{\mathcal{X}} I(x;\theta)\xi(x)dx.$$

The FIM establishes a connection between optimal experimental design and the Generalized Regression Model. The standard form of FIM under the normality hypothesis can be generalized to any probability distribution by including the EIM. By definition, the inverse of the FIM is asymptotically proportional to the variance and covariance matrix of estimators of *θ*, the parameters of the model. This matrix may depend on these parameters, so nominal values for them are necessary and therefore locally optimal designs can be obtained. By Carathéodory's theorem, it is known that for any design there is always another with the same information matrix of at most *m*(*m* + 1)/2 + 1 different points, where *m* is the number of unknown parameters to be estimated for the model *η*(*x*; *θ*) [17]. Therefore, it is sufficient to seek designs with finite support.

Optimization criteria express functions of the FIM that allow this matrix to be optimized in different ways. Consider the criterion function Φ as a real convex bounded function defined over the space of the information matrix <sup>M</sup> <sup>=</sup> {*M*(*ξ*) : *<sup>ξ</sup>* <sup>∈</sup> <sup>Ξ</sup>}. A design *<sup>ξ</sup>*<sup>∗</sup> will then be <sup>Φ</sup>-optimal if *<sup>ξ</sup>*<sup>∗</sup> <sup>=</sup> arg min*ξ*∈<sup>Ξ</sup> <sup>Φ</sup>(*M*(*ξ*; *<sup>θ</sup>*)). A number of studies, for example Chapter 10 of [18], give the criteria most commonly used in the literature. This paper uses the D-optimality criterion, whose goal is to minimize the volume of the confidence ellipsoid of *θ*ˆ, the estimators of *θ*. This criterion may be expressed by

$$\sup \Phi\_D(\mathcal{M}(\xi;\theta)) = \log |\mathcal{M}^{-1}(\xi;\theta)|.$$

In practice this criterion is equivalent to maximizing the determinant of the information matrix. The General Equivalence Theorem (see [19]) is a tool that allows optimality of a given design under a specific criterion to be checked. The sensitivity function *ϕ*(*x*; *ξ*, *θ*) is defined as a directional derivative

$$\varphi(\mathfrak{x};\mathfrak{f};\mathfrak{e}) = \lim\_{\mathfrak{a}\to 0^{+}} \frac{\partial}{\partial \mathfrak{a}} \Phi[\mathcal{M}((1-\mathfrak{a})\mathfrak{f} + \mathfrak{a}\mathfrak{f}\_{\mathfrak{a}};\mathfrak{e})]\_{\mathfrak{e}}$$

where ¯ *ξ<sup>x</sup>* is an arbitrary design centered on a point *x*. Given an optimal design, *ξ*∗, we find that *<sup>ϕ</sup>*(*x*; *<sup>ξ</sup>*∗, *<sup>θ</sup>*) <sup>≥</sup> 0, and the equality is found in the support points of the optimal design. The sensitivity function for the D-optimization criterion is given by

$$\varrho(\mathbf{x}; \underline{\mathfrak{z}}, \boldsymbol{\theta}) = m - \nu(\eta(\mathbf{x}; \boldsymbol{\theta})) f^T(\mathbf{x}; \boldsymbol{\theta}) \mathbf{M}^{-1}(\underline{\mathfrak{z}}; \boldsymbol{\theta}) f(\mathbf{x}; \boldsymbol{\theta}). \tag{3}$$

The efficiency allows any design *ξ* to be compared to the Φ-optimal design *ξ*∗,

$$\text{eff}\_{\Phi}(\boldsymbol{\xi}|\boldsymbol{\xi}^{\boldsymbol{\alpha}}) = \frac{\Phi(\boldsymbol{M}(\boldsymbol{\xi}^{\boldsymbol{\alpha}};\boldsymbol{\theta}))}{\Phi(\boldsymbol{M}(\boldsymbol{\xi};\boldsymbol{\theta}))}.$$

Also, if Φ is positively homogeneous, the value of the efficiency can be interpreted practically. If the efficiency value is 0.7, this means that the Φ-optimal design can be used to obtain the same information, or equivalently, the same statistical inference of the estimators of the model parameters, with a saving of 30% of the observations. For D-optimization criterion, which is positively homogeneous, D-efficiency is calculated as follows:

eff*D*(*ξ*|*ξ*∗) = <sup>|</sup>*M*(*ξ*; *<sup>θ</sup>*)<sup>|</sup> <sup>|</sup>*M*(*ξ*∗; *<sup>θ</sup>*)<sup>|</sup> 1/*m* . (4)

This expression will be termed "efficiency" from here on, as there is no possible confusion.

#### **3. Variance Structure and EIM for a Heteroscedastic Normal Distribution**

In most applications in the context of optimal experimental design, the homoscedastic normal distribution is used. However, when the response follows the gamma or the Poisson distribution the variance depends on the explanatory variable. To compare in a fair way with these distributions it is considered the heteroscedastic normal distribution with a variance structure given by

$$\text{Var}[y] = k \text{E}[y]^{2r},\tag{5}$$

where *<sup>k</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> and *<sup>r</sup>* <sup>∈</sup> <sup>R</sup> are constants and <sup>E</sup>[*y*] = *<sup>η</sup>*(*x*; *<sup>θ</sup>*). Thus, taking *<sup>k</sup>* <sup>=</sup> 1, for a value of *r* = 0.5 the variance structure for the heteroscedastic normal distribution is similar to that of the Poisson distribution (Var[*y*] = E[*y*]). On the other hand, with *k* = 1/*α* and *r* = 1, the structure of the variance for the heteroscedastic normal distribution is Var[*y*] = E[*y*] 2/*α*, similar to the variance of the gamma distribution, Γ(*α*, *β*), when parameter *α* is constant. Finally, the case *r* = 0 corresponds to the homoscedastic normal distribution.

Then, using (2), the EIM for the heteroscedastic normal distribution with variance given by (5) is

$$\nu\_N(\eta(\mathbf{x};\boldsymbol{\theta});r,k) = \frac{2r^2}{\eta(\mathbf{x};\boldsymbol{\theta})^2} + \frac{1}{k\eta(\mathbf{x};\boldsymbol{\theta})^{2r}}.$$

**Theorem 1.** *Let η*(*x*; *θ*) > 0 *be the function of some regression model, for any optimization criterion* Φ *based on the FIM, then the* Φ*-optimal designs for the heteroscedastic normal distribution with r* = 1 *in the variance defined in* (5) *and for the gamma distribution with constant α coincide. Also, the* Φ*-optimal design obtained is independent of α and k.*

**Proof.** Taking *r* = 1 in the variance given by (5), the EIM for the heteroscedastic normal distribution is *<sup>ν</sup>*<sup>N</sup> (*η*(*x*; *<sup>θ</sup>*)) = (2*<sup>k</sup>* <sup>+</sup> <sup>1</sup>)/(*kη*(*x*; *<sup>θ</sup>*)2), while the EIM for the gamma distribution is *ν*Γ(*η*(*x*; *θ*)) = *α*/*η*(*x*; *θ*)2, and so

$$\mathcal{M}\_{\mathcal{N}}(\xi;\theta) = \frac{2k+1}{k} \frac{1}{\eta(\mathbf{x};\theta)^2} f^T(\mathbf{x};\theta) f(\mathbf{x};\theta) \approx \frac{\mathfrak{a}}{\eta(\mathbf{x};\theta)^2} f^T(\mathbf{x};\theta) f(\mathbf{x};\theta) = \mathcal{M}\_{\Gamma}(\xi;\theta).$$

Therefore the Φ-optimal design calculated with any of the matrices will agree. Also, the parameters *k* and *α* are constants, multiplied in each expression of the FIM, and so do not affect Φ-optimal design.

The form of the EIMs of heteroscedastic normal (*r* = 0.5) and Poisson distribution are hardly proportional. Therefore, in this case, there is no possible similar result to Theorem 1.

#### **4. Linear Quadratic Model**

The linear quadratic model is considered in many studies which assume different probability distributions, such as gamma or Poisson distributions (for instance, refs. [1,4]). The regression function of the model is given by

$$\eta(\mathbf{x}; \boldsymbol{\theta}) = \theta\_0 + \theta\_1 \mathbf{x} + \theta\_2 \mathbf{x}^2, \ \mathbf{x} \in \mathcal{X}$$

The aim of this section is to provide D-optimal designs for this model when the response variable follows first a gamma and then a Poisson distribution. It also discusses the influence of misspecification for an assumed heteroscedastic normal distribution.

#### *4.1. Gamma Distribution*

Gamma models are suitable when the response is non-negative, continuous, skewed and heteroscedastic [7]. The introduction of the cited reference mentions several papers with real applications. From the point of view of optimal experimental design some papers could be cited, for example [6,20] for the case of multivariate gamma models, and [4] for the univariate case. In the present study, this last reference is revisited as an example of the applicability of the following results.

**Theorem 2.** *Let <sup>η</sup>*(*x*; *<sup>θ</sup>*) = *<sup>θ</sup>*<sup>0</sup> + *<sup>θ</sup>*1*<sup>x</sup>* + *<sup>θ</sup>*2*x*<sup>2</sup> + ... *<sup>θ</sup>px<sup>p</sup>* > <sup>0</sup> *be the function of a linear regression model of order <sup>p</sup>* <sup>≥</sup> <sup>1</sup>*, where <sup>x</sup> is defined on a design space* <sup>X</sup> = [*xl*, *xu*]*. If the response variable follows a gamma distribution with constant parameter α, the D-optimal design is supported in <sup>p</sup>* + <sup>1</sup> *equally weighted points with x*<sup>1</sup> = *xl and xp*+<sup>1</sup> = *xu. It can be expressed by*

$$\mathbf{f}\_{\Gamma}^{\*} = \begin{Bmatrix} \mathbf{x}\_{1} & \mathbf{x}\_{2} & \dots & \mathbf{x}\_{p} & \mathbf{x}\_{p+1} \\ 1/(p+1) & 1/(p+1) & \dots & 1/(p+1) & 1/(p+1) \end{Bmatrix} \text{.}$$

*For the linear quadratic model (p* = 2*), the D-optimal design is*

$$\mathbf{g}\_{\Gamma}^{\mathbf{x}\*} = \begin{Bmatrix} \mathbf{x}\_{l} & \mathbf{x}\_{2} & \mathbf{x}\_{u} \\ 1/3 & 1/3 & 1/3 \end{Bmatrix} \prime$$

*where x*<sup>2</sup> <sup>∈</sup> (*xl*, *xu*) *is a root of the linear quadratic equation*

$$(\theta\_1 + \theta\_2(\mathbf{x}\_l + \mathbf{x}\_u))\mathbf{x}\_2^2 - (2\theta\_2\mathbf{x}\_l\mathbf{x}\_u - 2\theta\_0)\mathbf{x}\_2 - (\theta\_0(\mathbf{x}\_l + \mathbf{x}\_u) + \theta\_1\mathbf{x}\_l\mathbf{x}\_u) = 0. \tag{6}$$

*Thus, it will be one of the solutions of*

$$\mathbf{x}\_{2} = \frac{2\theta\_{2}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{u}} - 2\theta\_{0} \pm \sqrt{(2\theta\_{2}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{u}} - 2\theta\_{0})^{2} + 4(\theta\_{1} + \theta\_{2}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{u}}))(\theta\_{0}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{u}}) + \theta\_{1}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{u}})}{2(\theta\_{1} + \theta\_{2}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{u}}))}.$$

**Proof.** Particularizing the sensitivity function given in (3) using the EIM for the gamma distribution (Table 1) gives

$$\varrho(\mathbf{x}; \underline{\mathfrak{z}}, \boldsymbol{\theta}) = (p+1) - \frac{\alpha}{\eta(\underline{\mathfrak{x}}; \boldsymbol{\theta})^2} f(\mathbf{x})^T \mathbf{M}^{-1}(\underline{\mathfrak{z}}; \boldsymbol{\theta}) f(\mathbf{x}).$$

By the General Equivalence Theorem, if *ξ*∗ <sup>Γ</sup> is the D-optimal design, *<sup>ϕ</sup>*(*x*; *<sup>ξ</sup>*<sup>∗</sup> <sup>Γ</sup>, *<sup>θ</sup>*) <sup>≥</sup> <sup>0</sup> must be satisfied for all *<sup>x</sup>* <sup>∈</sup> [*xl*, *xu*], and there must be equality in the support points of the design. It is, therefore, necessary to study the zeros of the function

$$\mathcal{g}(\mathbf{x}) = (p+1)\eta\left(\mathbf{x};\boldsymbol{\theta}\right)^2 - \alpha f(\mathbf{x})^T \mathbf{M}^{-1}(\boldsymbol{\xi}\_\Gamma^\*;\boldsymbol{\theta})f(\mathbf{x}),$$

which is a 2*p*-order polynomial and its zeros coincide with the zeros of *ϕ*(*x*; *ξ*<sup>∗</sup> <sup>Γ</sup>; *<sup>θ</sup>*). First, the number of support points of the D-optimal design must be greater or equal to the number of unknown parameters in the model, *m* = *p* + 1, in order for the FIM to be regular. Suppose, then, that the D-optimal design *ξ*∗ <sup>Γ</sup> has *<sup>p</sup>* <sup>+</sup> 2 support points. In this case, there will be at least *p* internal points with multiplicity two for the sensitivity function and its derivative to vanish, and the polynomial *g*(*x*) will have at least 2*p* + 2 roots, contradicting its order, which is 2*p*. Therefore, the D-optimal design cannot have more than nor fewer than *p* + 1 points, and so must have exactly *p* + 1 points. Now suppose that one extreme of X is not a support point of the design. Then it is assumed, without loss of generality, that the support points of the optimal design *<sup>x</sup>*1, ... , *xp*+<sup>1</sup> satisfies *xl* < *<sup>x</sup>*<sup>1</sup> < ... < *xp*+<sup>1</sup> = *xu*. The points *<sup>x</sup>*1, *<sup>x</sup>*2, ... , *xp* are roots of multiplicity 2 of *<sup>g</sup>*(*x*), and by Rolle's Theorem, there exist *<sup>c</sup>*<sup>1</sup> <sup>∈</sup> (*x*1, *<sup>x</sup>*2), *<sup>c</sup>*<sup>2</sup> <sup>∈</sup> (*x*2, *<sup>x</sup>*3), ... *cp* <sup>∈</sup> (*xp*, *xp*+1) such that *<sup>g</sup>* (*ci*) = 0, *<sup>i</sup>* = 1, ... *<sup>p</sup>*. Therefore, *g* (*x*) vanishes at 2*p* points, once again contradicting the order of the polynomial *g* (*x*), of order 2*<sup>p</sup>* <sup>−</sup> 1. By analogous reasoning, for the case *xl* <sup>=</sup> *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; ... <sup>&</sup>lt; *xp*+<sup>1</sup> <sup>&</sup>lt; *xu*, the conclusion is that the D-optimal design should have the two extremes in its support, and by the above, *p* − 1 internal points.

Finally, D-optimal design is equally weighted because the weights can be separated out in the optimization of the determinant in the way

$$|\mathcal{M}(\vec{\xi};\theta)| = \left(\prod\_{i=1}^{p+1} \nu(\eta(x\_i;\theta))\right) F(x\_1,\dots,x\_{p+1}) w\_1\dots w\_{p+1}$$

where

$$F(\mathbf{x}\_1, \dots, \mathbf{x}\_{p+1}) = \prod\_{\substack{i=1, j=2\\i$$

only depends on the support points. Thus, the maximum product of the *p* + 1 weights, which are restricted to being positive and summing to 1, is reached for *wi* = 1/(*<sup>p</sup>* + <sup>1</sup>).

For *<sup>p</sup>* = 2, the internal point of the design is found by solving, with *<sup>x</sup>*<sup>1</sup> = *xl* and *<sup>x</sup>*<sup>3</sup> = *xu*, the equation

$$\frac{\partial |M(\mathbf{j};\theta)|}{\partial \mathbf{x}\_2} = \frac{2(\mathbf{x}\_2 - \mathbf{x}\_l)(\mathbf{x}\_2 - \mathbf{x}\_u)(\mathbf{x}\_l - \mathbf{x}\_u)^2 a(\mathbf{x}\_2; \theta) w\_1 w\_2 w\_3}{\eta(\mathbf{x}\_2; \theta)^3 \eta(\mathbf{x}\_l; \theta)^2 \eta(\mathbf{x}\_u; \theta)^2} = 0 \tag{7}$$

where *<sup>a</sup>*(*x*2; *<sup>θ</sup>*)=(*θ*<sup>1</sup> + *<sup>θ</sup>*2(*xl* + *xu*))*x*<sup>2</sup> <sup>2</sup> <sup>−</sup> (2*θ*2*xlxu* <sup>−</sup>2*θ*0)*x*<sup>2</sup> <sup>−</sup> (*θ*0(*xl* <sup>+</sup> *xu*) + *<sup>θ</sup>*1*xlxu*). To solve Equation (7) is equivalent to solve *<sup>a</sup>*(*x*2) = 0, which is a linear quadratic equation with roots

$$\mathbf{x}\_{2} = \frac{2\theta\_{2}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{ul}} - 2\theta\_{0} \pm \sqrt{(2\theta\_{2}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{ul}} - 2\theta\_{0})^{2} + 4(\theta\_{1} + \theta\_{2}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{ul}}))(\theta\_{0}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{ul}}) + \theta\_{1}\mathbf{x}\_{l}\mathbf{x}\_{\mathrm{ul}})}{2(\theta\_{1} + \theta\_{2}(\mathbf{x}\_{l} + \mathbf{x}\_{\mathrm{ul}}))}.\text{}$$

By the previous results, only one of the two roots can be on the interval (*xl*, *xu*).

**Corollary 1.** *Let <sup>η</sup>*(*x*; *<sup>θ</sup>*) = *<sup>θ</sup>*<sup>0</sup> + *<sup>θ</sup>*1*<sup>x</sup>* + *<sup>θ</sup>*2*x*<sup>2</sup> + ... *<sup>θ</sup>px<sup>p</sup>* > <sup>0</sup> *be the function of a linear regression model of order <sup>p</sup>* <sup>≥</sup> <sup>1</sup>*, where <sup>x</sup> is defined on a design space* <sup>X</sup> = [*xl*, *xu*]*. If the response variable follows a heteroscedastic normal distribution, with r* = 1 *in the variance defined by Equation* (5)*, then ξ*∗ <sup>N</sup> <sup>=</sup> *<sup>ξ</sup>*<sup>∗</sup> Γ*.*

**Proof.** This is a direct consequence of Theorems 1 and 2.

**Corollary 2.** *By the hypothesis of Theorem 2, the following specific cases exist where the internal point of the design, x*2*, does not depend on the values of the parameters θ of the model:*

**Proof.** The cases can be computed by algebraic manipulation from Equation (6).

In [4], Bayesian, A- and D-optimal designs are computed for linear models assuming gamma distribution. In the case of linear quadratic model D-optimal designs are computed for different nominal values. Some of them are not equally weighted or even they are supported in two points (singular designs). This might seem in contradiction to Theorem 2 above. However, it happens only for the nominal values *<sup>θ</sup>*(0) for which the linear predictor *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 0 for, at least, one *<sup>x</sup>* ∈ X . If *<sup>η</sup>*(*x*; *<sup>θ</sup>*) = 0, a problem occurs in the definition of EIM (*ν*(*η*(*x*; *θ*)) = 1/*η*(*x*; *θ*)2). On the other hand, the case *η*(*x*; *θ*) < 0 does not make mathematical sense since *η*(*x*; *θ*)−<sup>1</sup> = *E*[*y*] = *α*/*β*, where *α*, *β* > 0 are the parameters of the gamma distribution (see Table 1).

For all nominal values *<sup>θ</sup>*(0) for which *<sup>η</sup>*(*x*; *<sup>θ</sup>*(0)) > 0, Theorem 2 can be applied to obtain D-optimal designs. Thus, both extremes of the design space are included in Doptimal designs, all of which are equally weighted, and the inner points, *x*2, are obtained by solving (6) (Table 2). The D-optimality condition is verified by the General Equivalence Theorem, through the sensitivity function (3). In addition, for the nominal values *<sup>θ</sup>*(0) = (0.3, <sup>−</sup>0.3, 0.3), the first condition of Corollary 2 is satisfied. Thus, it can be shown that the D-optimal design is supported in the midpoint of X , which agrees with the D-optimal design for a homoscedastic normal distribution.

**Table 2.** Locally D-optimal designs {*x*<sup>1</sup> <sup>=</sup> 0, *<sup>x</sup>*2, *<sup>x</sup>*<sup>3</sup> <sup>=</sup> <sup>1</sup>} equally weighted obtained for the linear quadratic regression model when the probability distribution of the response is gamma. The nominal values of the parameters of the model are those considered in [4].

*4.2. Poisson Distribution*

Generalized Linear Models for Poisson distribution are widely used in the literature. Special attention is paid to linear quadratic models in oncology [21–24]. A reference involving optimal designs and Poisson distribution is [1], where different linear regression models are considered.

**Theorem 3.** *Let <sup>η</sup>*(*x*; *<sup>θ</sup>*) = *<sup>θ</sup>*<sup>0</sup> + *<sup>θ</sup>*1*<sup>x</sup>* + *<sup>θ</sup>*2*x*<sup>2</sup> *be the function of the linear quadratic regression model, with <sup>x</sup> defined on the design space* <sup>X</sup> = [*xl*, *xu*]*, and the response variable following a Poisson distribution. Then, for the 3-point D-optimal design, we have the following sufficient conditions:*

*Also, if both extremes of* X *are included in the design, the internal point x*<sup>2</sup> *will be the solution, included in* X *, of the cubic equation*

$$\begin{aligned} -2\theta\_2 \mathbf{x}\_2^3 + [2\theta\_2(\mathbf{x}\_\mathsf{u} + \mathbf{x}\_l) - \theta\_1] \mathbf{x}\_2^2 + [\theta\_1(\mathbf{x}\_\mathsf{u} + \mathbf{x}\_l) - 2\mathbf{x}\_l \mathbf{x}\_\mathsf{u} \theta\_2 - 4] \mathbf{x}\_2 \\ + [2\mathbf{x}\_\mathsf{u} + \mathbf{x}\_l(2 - \mathbf{x}\_\mathsf{u}\theta\_1)] = 0. \end{aligned} \tag{8}$$

**Proof.** Consider the 3-point D-optimal design

$$
\xi\_{\mathcal{P}}^{\mathbf{x}\*} = \begin{Bmatrix} \mathbf{x}\_1 & \mathbf{x}\_2 & \mathbf{x}\_3 \\ 1/3 & 1/3 & 1/3 \end{Bmatrix},
$$

with *xl* <sup>≤</sup> *<sup>x</sup>*<sup>1</sup> <sup>&</sup>lt; *<sup>x</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>x</sup>*<sup>3</sup> <sup>≤</sup> *xu*. The design is equally weighted because the weights can be separated out in the optimization of the determinant (see Proof of Theorem 2).

The explicit expression of the derivative with respect to *x*<sup>1</sup> is

$$\begin{split} \frac{\partial |\mathcal{M}(\frac{\chi}{6};\theta)|}{\partial \mathbf{x}\_{1}} &= \frac{1}{27} \exp\left\{ \sum\_{i=1}^{3} \eta \left( \mathbf{x}\_{i};\theta \right) \right\} \left( \mathbf{x}\_{2} - \mathbf{x}\_{1} \right) \left( \mathbf{x}\_{3} - \mathbf{x}\_{1} \right) \left( \mathbf{x}\_{3} - \mathbf{x}\_{2} \right)^{2} \\ &\times \left( \left( 4\mathbf{x}\_{1} - 2\mathbf{x}\_{2} - 2\mathbf{x}\_{3} \right) + \left( \mathbf{x}\_{2} - \mathbf{x}\_{1} \right) \left( \mathbf{x}\_{3} - \mathbf{x}\_{1} \right) \left( \theta\_{1} + 2\mathbf{x}\_{1} \theta\_{2} \right) \right), \end{split}$$

If *<sup>∂</sup>*|*M*(*ξ*; *<sup>θ</sup>*)|/*∂x*<sup>1</sup> <sup>&</sup>lt; 0 on [*xl*, *<sup>x</sup>*2), then the maximum of the determinant will be reached at *<sup>x</sup>*<sup>1</sup> = *xl*. Thus

$$\frac{\partial |\mathcal{M}(\vec{\xi};\theta)|}{\partial \mathbf{x}\_{1}} < 0 \iff \theta\_{1} + 2\mathbf{x}\_{1}\theta\_{2} < \frac{2\mathbf{x}\_{2} + 2\mathbf{x}\_{3} - 4\mathbf{x}\_{1}}{(\mathbf{x}\_{2} - \mathbf{x}\_{1})(\mathbf{x}\_{3} - \mathbf{x}\_{1})}.\tag{9}$$

If we consider *<sup>θ</sup>*<sup>2</sup> <sup>&</sup>lt; 0, we have *<sup>θ</sup>*<sup>1</sup> <sup>+</sup> <sup>2</sup>*x*1*θ*<sup>2</sup> <sup>≤</sup> *<sup>θ</sup>*<sup>1</sup> <sup>+</sup> <sup>2</sup>*xlθ*<sup>2</sup> and the inequalities

$$\frac{2\mathbf{x}\_2 + 2\mathbf{x}\_3 - 4\mathbf{x}\_1}{(\mathbf{x}\_2 - \mathbf{x}\_1)(\mathbf{x}\_3 - \mathbf{x}\_1)} > \frac{4(\mathbf{x}\_2 - \mathbf{x}\_1)}{(\mathbf{x}\_2 - \mathbf{x}\_1)(\mathbf{x}\_3 - \mathbf{x}\_1)} = \frac{4}{(\mathbf{x}\_3 - \mathbf{x}\_1)} > \frac{4}{(\mathbf{x}\_\mu - \mathbf{x}\_1)}.$$

are satisfied.

Therefore, the inequality (9) is true if the following is satisfied

$$
\theta\_1 + 2\mathbf{x}\_l \theta\_2 < \frac{4}{\left(\mathbf{x}\_u - \mathbf{x}\_l\right)^2}.
$$

Also,

$$\begin{split} \frac{\partial |M(\frac{\eta}{3};\theta)|}{\partial \mathbf{x}\_{3}} &= \frac{1}{27} \exp\left\{ \sum\_{i=1}^{3} \eta \left( \mathbf{x}\_{i};\theta \right) \right\} (\mathbf{x}\_{3} - \mathbf{x}\_{1}) (\mathbf{x}\_{3} - \mathbf{x}\_{2}) (\mathbf{x}\_{2} - \mathbf{x}\_{1})^{2} \\ &\times \left( (4\mathbf{x}\_{3} - 2\mathbf{x}\_{1} - 2\mathbf{x}\_{2}) + (\mathbf{x}\_{3} - \mathbf{x}\_{1})(\mathbf{x}\_{3} - \mathbf{x}\_{2}) (\theta\_{1} + 2\mathbf{x}\_{3}\theta\_{2}) \right), \end{split}$$

and if *<sup>∂</sup>*|*M*(*ξ*; *<sup>θ</sup>*)|/*∂x*<sup>3</sup> <sup>&</sup>gt; 0 on (*x*2, *xu*] the maximum will be found at *<sup>x</sup>*<sup>3</sup> <sup>=</sup> *xu*. This gives

$$\frac{\partial |\mathcal{M}(\mathbf{\hat{y}};\theta)|}{\partial \mathbf{x}\_3} > 0 \iff \frac{2\mathbf{x}\_1 + 2\mathbf{x}\_2 - 4\mathbf{x}\_3}{(\mathbf{x}\_3 - \mathbf{x}\_1)(\mathbf{x}\_3 - \mathbf{x}\_2)} < \theta\_1 + 2\mathbf{x}\_3\theta\_2. \tag{10}$$

If *<sup>θ</sup>*<sup>2</sup> <sup>&</sup>lt; 0, we have *<sup>θ</sup>*<sup>1</sup> <sup>+</sup> <sup>2</sup>*xuθ*<sup>2</sup> <sup>≤</sup> *<sup>θ</sup>*<sup>1</sup> <sup>+</sup> <sup>2</sup>*x*3*θ*2. Thus,

$$\frac{2x\_1 + 2x\_2 - 4x\_3}{(x\_3 - x\_1)(x\_3 - x\_2)} < 0\_{\prime\prime}$$

and so if 0 < *<sup>θ</sup>*<sup>1</sup> + <sup>2</sup>*xuθ*2, the inequality in (10) is satisfied.

Finally, if *xl* and *xu* are in the support of the design, like in Theorem 2, the internal point will be a solution of the equation *<sup>∂</sup>*|*M*(*ξ*; *<sup>θ</sup>*)|/*∂x*<sup>2</sup> <sup>=</sup> 0, which is equivalent to the cubic equation given by (8).

As mentioned above, the linear quadratic model plays an important role in oncology, and optimal experimental design has an important practical role in determining the best doses for carrying out the experiment and fitting the model. To illustrate the previous result, we consider the example in [1] where the response variable *y* explains the number of living cells in a system and the explanatory variable *x* is the dose of an injected oncology drug. Hence, the expected number of living cells for any dose *xi* is given by

$$
\lambda\_i = E[y\_i] = e^{\theta\_0 + \theta\_1 x\_i + \theta\_2 x\_i^2} \quad \text{x}\_i \ge 0.1
$$

From the context of the problem, the relationship between *x* and *y* must be inverse: the higher the dose inoculated the lower the number of living cells and vice versa. For the examples, *<sup>θ</sup>*<sup>1</sup> <sup>≤</sup> 0 and *<sup>θ</sup>*<sup>2</sup> <sup>&</sup>lt; 0 are considered to satisfy this relationship for all *<sup>x</sup>* ∈ X . Furthermore, to consider a high dose would not be realistic, as the number of living cells could be very low and might compromise the survival of the system. Let *λ<sup>c</sup>* = *eθ*<sup>0</sup> be the mean of the number of surviving cells for a control dose (*x* = 0). Then, the expected survival proportion for any dose *xi* is *<sup>λ</sup>i*/*λ<sup>c</sup>* <sup>≥</sup> *<sup>c</sup>*, where *<sup>c</sup>* <sup>∈</sup> (0, 1] is the minimal survival proportion. The value of *c* is a characteristic for each system and for the context of the problem. For this study we consider *c* = 0.4. When *θ*<sup>2</sup> <sup>1</sup>/*θ*<sup>2</sup> ≥ −4 log *c*, the survival proportion is not less than the minimal survival proportion in the design space <sup>X</sup> = [0, *xu*], where *xu* is expressed as a function of the parameters of the model (see details in [1]).

Based on the above, the first condition of the Theorem 3 is satisfied and therefore a control dose *x* = 0 is always included on the D-optimal design. Table 3 shows D-optimal designs when the response variable follows heteroscedastic normal (with *r* = 0.5 in (5)) or Poisson distributions. The nominal values considered fulfill the relationship in [1]. Moreover, all D-optimal designs are supported on the upper extreme, so only the inner points *x*<sup>2</sup> of D-optimal designs are shown. For a Poisson distribution the point *x*<sup>2</sup> may be computed by solving Equation (8) of Theorem 3. Finally, an efficiency study is carried out. The efficiencies of the designs are calculated by adapting (4) as

$$\text{eff}\_D(\ulcorner\_A | \ulcorner\_T) = \left(\frac{|\mathcal{M}\_T(\ulcorner\_A; \theta)|}{|\mathcal{M}\_T(\ulcorner\_T; \theta)|}\right)^{1/m} \urcorner \tag{11}$$

where *ξ<sup>A</sup>* is the D-optimal design for the probability distribution assumed by the researcher (for this example, heteroscedastic normal with *<sup>r</sup>* = 0.5), while *<sup>ξ</sup><sup>T</sup>* is the D-optimal design and *M<sup>T</sup>* is the FIM, both for the true probability distribution (in this example, a Poisson distribution). The last column of Table 3 shows that efficiencies. Unlike the results obtained for the gamma distribution, where the D-optimal designs coincide with the heteroscedastic normal distribution when the relationship between mean and variance agrees (Corollary 1), there is a non-negligible loss efficiency, around 20% or more, in this case. It is noteworthy that the inner point of the Poisson distribution is lower than that for a heteroscedastic normal distribution for the designs computed.

**Table 3.** Locally D-optimal designs {*x*<sup>1</sup> <sup>=</sup> 0, *<sup>x</sup>*2, *<sup>x</sup>*<sup>3</sup> <sup>=</sup> *xu*} equally weighted obtained for the linear quadratic model when the probability distribution of the response is heteroscedastic normal or Poisson. The nominal values are *<sup>θ</sup>*(0) = (0.95, <sup>−</sup>1, *<sup>θ</sup>*2) and the minimal survival fraction is *<sup>c</sup>* <sup>=</sup> 0.4. The last column shows the efficiency when comparing the design for the heteroscedastic normal distribution with *r* = 0.5 to the Poisson distribution.

#### **5. Extended Hill Model**

The Hill model is a dose-response model commonly used in practice to describe the relationship between the concentration of a drug and its effect. Several papers [12,14,15] have addressed this issue from the point of view of optimal experimental design. This model may explain both discrete and continuous responses, such as counting cells [25] or the effect of a drug on cell growth [13], among many others. Here we focus on the 4-parameter Hill model.

If we consider *x* to be the dose of an administered drug, the function of the regression model which explains the effect can be expressed as

$$\eta(\mathbf{x}; \boldsymbol{\theta}) = \frac{(E\_{\rm con} - b) \left(\frac{\mathbf{x}}{\|\mathbf{C}\_{\rm S0}\|}\right)^{s}}{1 + \left(\frac{\mathbf{x}}{\|\mathbf{C}\_{\rm S0}\|}\right)^{s}} + b,\tag{12}$$

where *<sup>θ</sup>* = (*Econ*, *<sup>b</sup>*, *IC*50,*s*) are the parameters to be estimated. The parameter *Econ* is the effect on the control, i.e., where there is no dose. The parameter *b* corresponds to the asymptotic value of the response when the concentration of the drug tends to infinity and *IC*<sup>50</sup> corresponds to the dose at which a response would be found equal to the middle of the effect range, *Econ* <sup>−</sup> *<sup>b</sup>*. Finally, the parameter *<sup>s</sup>* is a form parameter: if *<sup>s</sup>* <sup>&</sup>gt; 0, *<sup>η</sup>*(*x*; *<sup>θ</sup>*) will be strictly increasing, and if *s* < 0, strictly decreasing. Thus, when the parameter *b* > 0 and *s* < 0 the drug has an inhibitory effect where *b* implies that the whole cell population is not destroyed, as shown in Figure 1. This is the case considered in this paper. Here it is studied from two perspectives simultaneously, where the gamma, or Poisson, is the true distribution of the response variable and the practitioner assumes a heteroscedastic normal distribution with the variance structure given by (5).

**Figure 1.** Graph of the regression function *η*(*x*; *θ*) of the 4-parameter Hill model with nominal values corresponding to the drug TMTX, shown in Table 4.

Ref. [13] bring together different maximum likelihood estimations of the parameters of (12) for different types of drugs. Table 4 shows these nominal values and the 4-point D-optimal designs obtained for different probability distributions of the response variable: *ξ*<sup>Γ</sup> (gamma distribution with constant *α*), *ξ*<sup>N</sup> (heteroscedastic normal distribution with variance structure given by (5)) and *<sup>ξ</sup>*<sup>P</sup> (Poisson distribution). By Theorem 1, when *<sup>r</sup>* <sup>=</sup> <sup>1</sup> in (5) *<sup>ξ</sup>*<sup>Γ</sup> and *<sup>ξ</sup>*<sup>N</sup> coincide. However, the designs *<sup>ξ</sup>*<sup>P</sup> and *<sup>ξ</sup>*<sup>N</sup> with *<sup>r</sup>* <sup>=</sup> 0.5 are distinct, even though both comparisons show a similar relationship between the mean and the variance. Table 4 shows only the inner points (intermediate doses) of the D-optimal designs, as the extremes of the design space <sup>X</sup> = [0, *Dmax*] are included in all the cases studied. The maximum dosage *Dmax* was given by the value 1000 · *IC*50, except for the drug AG2009, since the authors considered this dosage to be impractical. It can be seen that the D-optimal design leads to experimenting with three very low doses, and at the maximum dose (*Dmax*), except for drug AG2009, where the doses are more spread out.

**Table 4.** Locally D-optimal designs {*x*<sup>1</sup> <sup>=</sup> 0, *<sup>x</sup>*2, *<sup>x</sup>*3, *<sup>x</sup>*<sup>4</sup> <sup>=</sup> *Dmax*} equally weighted for the 4-parameter Hill model for different drugs and probability distributions. The nominal values *<sup>E</sup>*(0) *con* <sup>=</sup> 1.70 and *<sup>b</sup>*(0) <sup>=</sup> 0.137 were considered for all the drugs. Columns 2–4 show the nominal values of the parameters and columns 5–10 show the internal points of the D-optimal designs *ξ*<sup>N</sup> , *ξ*<sup>Γ</sup> and *ξ*P. The last column shows the efficiency when comparing the design for the heteroscedastic normal distribution with *r* = 0.5 to the Poisson distribution.

The last column of Table 4 shows that the efficiency computed by (11) of the D-optimal designs when a heteroscedastic normal distribution with *r* = 0.5 is assumed with respect to the Poisson distribution, is around 73%, except for the drug AG2009, whose efficiency is higher. Again, in this practical case there is a considerable loss of efficiency in estimating the model parameters, with regard to misspecification in the probability distribution. All D-optimal designs in this section have been computed using the Wynn-Fedorov's algorithm [26].

#### *Sensitivity Analysis*

The main aim of this section is to study the effect of the relationship between E[*y*] and Var[*y*], characterized by the parameter *r* in (5), on the efficiency. So, a sensitivity analysis of this parameter is done. Ref. [11] studies the influence of misspecification in the structure of the variance in an analysis carried out for the gamma distribution and the heteroscedastic normal distribution separately. Here, a similar study was carried out with a point of view in which a practitioner considers a heteroscedastic normal response, but the true distribution of the response is gamma, or Poisson. For both distributions, efficiencies, using (11), are computed by comparing D-optimal designs for heteroscedastic normal distribution with the D-optimal design for the true probability distribution, as a function of the values of *r*.

The efficiencies achieved for different drugs are shown in Figure 2. It can be seen how, when the true distribution is gamma (Figure 2a), the efficiency is 1 for *r* = 1 (dot), since in this case the designs coincide as proven in Theorem 1. However, when the true distribution is the Poisson (Figure 2b), maximum efficiency is not obtained for *r* = 0.5 (dots), as might be expected. It is achieved in this case for negative values of *r*, close to *r* = 0, and so it would have been better, in terms of efficiency, for the practitioner to have assumed the homoscedastic normal distribution rather than heteroscedastic normal with *r* = 0.5. Furthermore, it does not reach the value 1 for any value of *r*. Finally, for this model in the neighborhood of *r* = 0 (homoscedastic normal distribution) opposite effects are produced on the efficiency for each of the distributions: greater efficiencies when the true distribution is the Poisson, and lower in the case of the gamma. It is important to highlight that there is no analytic explanation for this effect, and it is motivated by the model and nominal values.

**Figure 2.** Efficiencies when comparing the designs obtained for the heteroscedastic normal distribution with variance given by (5) and different values of *r* when the true distributions are the gamma (**a**) and Poisson (**b**), for the 4-parameter Hill model using the nominal values of Table 4. The graphs for the drugs MTX, AG2032 and AG2034 are similar to the graph for the drug TMTX.

For that, the effect of *r* on the trend of the efficiency is studied depending on the values taken by *η*(*x*; *θ*). Again, in this analysis, it is assumed that *y* follows a heteroscedastic normal distribution with variance structure given by (5) when the distribution of *y* is the Poisson or the gamma and a misspecification takes place.

First, for sufficiently large values of*<sup>r</sup>* and *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>&</sup>gt; 1, <sup>∀</sup>*<sup>x</sup>* ∈ X , given that *<sup>η</sup>*(*x*; *<sup>θ</sup>*)−2*r*/*<sup>k</sup>* <sup>≈</sup> 0, we have

$$\begin{aligned} \nu\_N(\eta(\mathbf{x};\boldsymbol{\theta}),r,k) &= 2r^2/\eta(\mathbf{x};\boldsymbol{\theta})^2 + \eta(\mathbf{x};\boldsymbol{\theta})^{-2r}/k \approx 2r^2/\eta(\mathbf{x};\boldsymbol{\theta})^2 \\ &\propto \boldsymbol{\alpha}/\eta(\mathbf{x};\boldsymbol{\theta})^2 = \nu\_\Gamma(\eta(\mathbf{x};\boldsymbol{\theta})). \end{aligned}$$

Thus, when the true probability distribution is a gamma distribution, Figure 3a (solid line) shows how, on increasing the value of *r* the efficiency tends to 1. On the other hand, the lower the value of *<sup>r</sup>*, the greater the difference between *<sup>ν</sup>*<sup>N</sup> (*η*(*x*; *<sup>θ</sup>*);*r*, *<sup>k</sup>*) and *ν*Γ(*η*(*x*; *θ*)), therefore the efficiency tends to 0 as can be seen in Figure 3a. However, if 0 <sup>&</sup>lt; *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 1 (dashed line), <sup>∀</sup>*<sup>x</sup>* ∈ X , the effect of *<sup>r</sup>* on the trend of the efficiencies of the designs obtained for the heteroscedastic normal distribution when the true distribution is a gamma distribution is the opposite. As it is shown in Figure 3a, if *r* increases, the efficiency tends to 0, and if *r* decreases, the efficiency tends to 1.

**Figure 3.** Study of efficiency trend when comparing the designs obtained for the heteroscedastic normal distribution as a function of the parameter *r* of (5) considering the gamma (**a**) or Poisson (**b**) distributions as the true distributions, for the 4-parameter Hill model in <sup>X</sup> = [0, 1500]. Solid lines assume *Econ* = 10 and *b* = 1 (therefore *η*(*x*; *θ*) > 1), while for dashed lines *Econ* = 1 and *b* = 0.1 (0 <sup>&</sup>lt; *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 1). Both cases use the nominal values *IC*<sup>50</sup> <sup>=</sup> 550 and *<sup>s</sup>* <sup>=</sup> <sup>−</sup>2.

When the true distribution is a Poisson distribution there is no direct comparison between its EIM and the EIM of a heteroscedastic normal distribution. However, it can be seen in Figure 3b how the efficiency reaches a maximum for a particular value of *r* and loses efficiency for values away from that value. This is because the study looks at values of *<sup>s</sup>* <sup>&</sup>lt; 0, and so *<sup>η</sup>*(*x*; *<sup>θ</sup>*) and *<sup>ν</sup>*P(*η*(*x*; *<sup>θ</sup>*)) are monotonic. Therefore, the maximum efficiency is at the value of *<sup>r</sup>* where the distance between *<sup>ν</sup>*<sup>N</sup> (*η*(*x*; *<sup>θ</sup>*)) and *<sup>ν</sup>*P(*η*(*x*; *<sup>θ</sup>*)) is minimal (independently of whether *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>&</sup>gt; 1 or 0 <sup>&</sup>lt; *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 1). Although the 4-parameter Hill model is taken as an example and the graphs in Figure 3 are obtained based on that model, the whole study on the trend of *r* on the efficiency is general for any regression function satisfying the inequalities.

Finally, it is interesting to point out differences between the graphs in Figures 2 and 3. First, the trends in the efficiencies as a function of *r* do not coincide. This is because, for the drugs in the study for the 4-parameter Hill model, the inequalities *η*(*x*; *θ*) > 1 or <sup>0</sup> <sup>&</sup>lt; *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 1 are not satisfied on the design spaces considered in the examples. Secondly, when the true distribution is the Poisson distribution, maximum efficiency in Figure 3b is obtained for a value close to *<sup>r</sup>* <sup>=</sup> 0 when 0 <sup>&</sup>lt; *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>≤</sup> 1, as also it takes place in Figure 2b, while for the case with *<sup>η</sup>*(*x*; *<sup>θ</sup>*) <sup>&</sup>gt; 1 maximum efficiency is obtained close to *<sup>r</sup>* <sup>=</sup> <sup>−</sup>3, i.e., for the same model, the nominal values defined by the context of the problem affect the loss of efficiency.

#### **6. Summary and Conclusions**

This study has been carried out to analyze the effect of misspecification in the probability distribution of the response variable. We measure that effect by calculating the efficiency of the optimal design obtained with an assumed or working distribution compared to that obtained with the true probability distribution. The typical case is when a researcher assumes a normal distribution, even a heteroscedastic one, for the response variable of his or her problem, but at a greater depth, another distribution is more appropriate, for example a gamma (or Poisson) distribution. When there is misspecification in the probability distribution, there is a loss of efficiency which depends both on the assumed probability distribution and on the regression function *η*(*x*; *θ*).

We provide some theoretical results, as well as practical ones. The first is quite general, valid for any regression function and any criterion based on FIM which guarantees that there is no loss of efficiency when the response variable follows a gamma distribution, and there is assumed to be a heteroscedastic normal distribution with *r* = 1 in the variance structure given by (5). For the linear quadratic model, analytical results are obtained on computing the optimal design for Poisson and gamma distributions. These theoretical results have been used in real applications from the literature, providing designs useful for practitioners.

Finally, the 4-parameter Hill model was used to illustrate and quantify the loss of efficiency. Assuming a heteroscedastic normal distribution, taking values close to *r* = 0 in (5), between about 18% and 25% efficiency is lost for all the drugs looked at in the study when the true distribution is a gamma distribution. Thus, in this case, the usual assumption of normality and homoscedasticity (*r* = 0) of the response variable is not a good option. However, when the true distribution is the Poisson, the loss of efficiency is less severe, reaching maximum values of efficiency for values close to *r* = 0 for all the drugs. This is a striking case, as one might expect maximum efficiency to be achieved at the value *r* = 0.5, which leads to the same relationship between the mean and the variance for the heteroscedastic normal and the Poisson distributions.

It is worth finishing this paper by mentioning that the EIM is an essential tool, as it collects information both about the regression function and the probability distribution of the response variable. As already mentioned, to assume the homoscedastic normal distribution when obtaining optimal designs may lead to a great loss of efficiency. Nonetheless, the examples given show that this will depend on the true distribution of the response variable and on the model function chosen. The existence of uncertainty about the probability distribution of the response variable will therefore lead to the future goal of obtaining robust designs to reduce this uncertainty.

**Author Contributions:** Conceptualization, M.A.-S., V.C.-A. and S.P.-C.; methodology, M.A.-S., V.C.- A. and S.P.-C.; software, M.A.-S., V.C.-A. and S.P.-C.; formal analysis, M.A.-S., V.C.-A. and S.P.-C.; investigation, M.A.-S., V.C.-A. and S.P.-C.; writing—original draft preparation, M.A.-S., V.C.-A. and S.P.-C.; writing—review and editing, M.A.-S., V.C.-A. and S.P.-C. All authors have read and agreed to the published version of the manuscript.

**Funding:** All authors were sponsored by Ministerio de Economía y Competitividad and fondos FEDER MTM2016-80539-C2-1-R and by Junta de Comunidades de Castilla-La Mancha SBPLY/17/ 180501/000380.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **Optimal Experimental Design for Parametric Identification of the Electrical Behaviour of Bioelectrodes and Biological Tissues**

**Àngela Sebastià Bargues 1, José-Luis Polo Sanz <sup>2</sup> and Raúl Martín Martín 1,\***

**Abstract:** The electrical behaviour of a system, such as an electrode–tissue interface (ETI) or a biological tissue, can be used for its characterization. One way of accomplishing this goal consists of measuring the electrical impedance, that is, the opposition that a system exhibits to an alternating current flow as a function of frequency. Subsequently, experimental impedance data are fitted to an electrical equivalent circuit (EEC model) whose parameters can be correlated with the electrode processes occurring in the ETI or with the physiological state of a tissue. The EEC used in this paper is a reasonable approach for simple bio-electrodes or cell membranes, assuming ideal capacitances. We use the theory of optimal experimental design to identify the frequencies in which the impedance is measured, as well as the number of measurement repetitions, in such a way that the EEC parameters can be optimally estimated. Specifically, we calculate approximate and exact D-optimal designs by optimizing the determinant of the information matrix by adapting two of the most algorithms that are routinely used nowadays (REX random exchange algorithm and KL exchange algorithm). The D-efficiency of the optimal designs provided by the algorithms was compared with the design commonly used by experimenters and it is shown that the precision of the parameter estimates can be increased.

**Keywords:** optimal experimental design; bioimpedance; impedance spectroscopy; algorithm

#### **1. Introduction**

The electrical behaviour of cells, biological tissues, or electrode-tissue interfaces can be used for their characterization [1–4]. One way to do this is by measuring the electrical impedance, that is, the opposition that the system under study exhibits to an alternating current flow as a function of frequency. Subsequently, experimental impedance data are fitted to an electrical equivalent circuit (EEC), which is a mathematical model that approximates the electrical behaviour of the system under study. Thus, the EEC parameters give information about the physiological state of a tissue or about the electrode processes occurring in the electrode-tissue interface (ETI). For instance, it is useful to know the electrical/electrochemical behaviour of the ETI plays a significant role in the biopotential measurements or also in the propagation of an applied stimuli. It should be mentioned that the experimental impedance data are analysed by using an impedance function. This impedance function can be proposed from a plausible physical theory (that predicts the impedance) or from an EEC, that aids in the visualization of the physical processes occurring in the system under study [1].

In this paper, we work with linear time-invariant (LTI) circuits, that is, LTI systems that fulfil the properties of linearity and invariance in time, meaning that the superposition and proportionality principles are hold. From the study of the circuit theory, a LTI circuit in which the input *i*(*t*) is an electrical current (amperes), the output *v*(*t*) is a voltage (volts),

**Citation:** Sebastià Bargues, À.; Polo Sanz, J.-L.; Martín Martín, R. Optimal Experimental Design for Parametric Identification of the Electrical Behaviour of Bioelectrodes and Biological Tissues. *Mathematics* **2022**, *10*, 837. https://doi.org/10.3390/ math10050837

Academic Editor: Christophe Guyeux

Received: 30 January 2022 Accepted: 2 March 2022 Published: 6 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**<sup>\*</sup>** Correspondence<EMAIL_ADDRESS>
and *t* is the time (seconds) can be described by an ordinary, linear differential equation (ODE) with constant coefficients *ai*, *bj* <sup>∈</sup> <sup>R</sup> and order *<sup>n</sup>* for *<sup>n</sup>*, *<sup>m</sup>* <sup>∈</sup> <sup>N</sup> as [5]:

$$\frac{d^n v(t)}{dt^n} + a\_{n-1} \frac{d^{n-1} v(t)}{dt^{n-1}} + \dots + a\_0 v(t) = b\_m \frac{d^m i(t)}{dt^m} + b\_{m-1} \frac{d^{m-1} i(t)}{dt^{m-1}} + \dots + b\_0 i(t) \tag{1}$$

Without loss of generality, we have assumed that *an* = 1.

Let us now consider the input *i*(*t*) = *est*, where *s* is a time-independent parameter. It can be shown that the output (refer to Equation (1)) can be written as *v*(*t*) = *Z*(*s*)*est*, where *Z*(*s*) is an impedance function (in units of ohms, Ω), which relates the voltage to the current:

$$Z(s) = \frac{b\_{m}s^{m} + b\_{m-1}s^{m-1} + \dots + b\_{0}}{s^{n} + a\_{n-1}s^{n-1} + \dots + a\_{0}} \tag{2}$$

The input given by *est* is extremely versatile. It allows us to obtain the forced response to a constant input (*s* = 0) or to a sinusoidal input (*s* = *jω*, where *j* is the imaginary unit and *ω* is the angular frequency in rad/s with *ω* = 2*π f* , where *f* is frequency in Hz). This latter response is known as frequency response of systems, that is, the response of a system to input sinusoids of different frequencies. Moreover, *est* constitutes a basis function involving the Fourier series and the Fourier and Laplace transforms [5].

Note in particular that the parameter *s* in (2) can be interpreted as the variable that appears in the Laplace transform or as the differential operator ( *<sup>d</sup> dt* , typically denoted by *p*) [6]. When *s* = *jω*, the circuit is operating at sinusoidal steady-state. Specifically, by considering the Laplace transform of Equation (1), assuming zero initial conditions, we can write *V*(*s*) = *Z*(*s*)*I*(*s*), where *I*(*s*) and *V*(*s*) are the Laplace transforms of the input *i*(*t*) and the output *v*(*t*), respectively. Interestingly, the unit-impulse (*i*(*t*) = *δ*(*t*), *I*(*s*) = 1, where *δ*(*t*) is Dirac delta function) response is *V*(*s*) = *Z*(*s*), that is, the Laplace transfom of the voltage output equal the impedance *Z*(*s*).

Electrical impedance has been usually defined in the context of a single sinusoidal signal and phasor analysis [6]. The impedance *Z*(*jω*) at the frequency *ω*, is a complex number whose magnitude and phase are <sup>|</sup>*Z*(*jω*)<sup>|</sup> and arg[*Z*(*jω*)], respectively. Equivalently, *<sup>Z</sup>*(*jω*) can also be expressed in terms of its real and imaginary components, that is, the resistance *R*(*ω*) and the reactance *X*(*ω*), respectively.

$$Z(j\omega) = |Z(j\omega)|\_{\text{arg}[Z(j\omega)]} = R(\omega) + jX(\omega) \tag{3}$$

Note that for physical frequencies (*s* = *jω*), the impedance in Equation (2) can be written according to Equation (3). Hereinafter, the impedance *Z*(*jω*) is also denoted by *Z*(*ω*).

Regarding the sinusoidal signal, it should be mentioned that the square of the rms (root mean square) value of a sinusoidal signal is the average power associated with it, that is, the average power is concentrated at frequency *ω*. Importantly, the average power associated with a periodic signal is the sum of the squared rms values (sum of the average powers) of all its components (Parseval's theorem for periodic signals). Therefore, the power spectral density (which, when integrated over the whole spectrum frequency, gives the total average power) for a periodic signal (power signal) consists of a series of impulses [7].

In this paper, we focus on finding the input frequencies to the circuit which maximises the amount of information obtained from an experiment, given that the true system is a priori known. Optimal input designs concern finding the input signal, which assures that the estimates become as good as possible. The use of optimal input signals will increase the precision of the parameter estimators [8]. This problem, called optimal input design, arises from the works of Mehra [9,10]. The theory of optimal design of experiments for the construction of optimal input signals in control theory, involving the frequency domain, can be found in [11–18]. However, its use for electrical impedance models has been scarcely reported in the scientific literature. Considerations of optimal multisine input signals are analysed in [19–21]. Specifically, a D-optimal multisine excitation for broadband impedance spectroscopy measurements is proposed in [20]. These optimal designs are based on simple first-order LTI circuits, that is, they are described by an ODE of order 1. The objective of this work is to design the D-optimal frequencies in which to carry out the electrical impedance measurements to achieve the best statistical estimates. Specifically, we calculate approximate and exact optimal designs by optimizing the determinant of the information matrix by adapting the REX random exchange algorithm and KL exchange algorithm. To the best of the authors' knowledge, there is no previous report of modifying those algorithms to calculate D-optimal designs in the frequency domain. In Section 2, we provide an introduction to the D-optimal input design applied to the study of a simple impedance model, which approaches the electrical behaviour of basic bio-electrodes or cell membranes. The definition of the Fisher Information Matrix (FIM) and the spectral density function are presented. Section 3 is devoted to the D-optimization criterion, the directional derivative, the equivalence theorem and the two algorithms adapted. Finally, in Section 4, we include a real application. The results obtained from an experimental test carried out with the two algorithms adapted to compare the efficiency of the design commonly used by experimenters with the D-optimal designs obtained from the algorithms are also discussed in this section.

#### **2. Optimal Experimental Design for the Identification of LTI Systems**

In general, the procedure for the identification of a system can be divided into three steps. The first one will be the experimental design, that is, planning and selecting the input signals to be introduced into the model. The set where to take the observations is known as the design space and is denoted by *χ*. From now on we will assume that it is a compact set. Once the input signals have been selected within the design space, the second step will be to collect data. The third step will be to estimate the parameters of the model.

#### *2.1. Design of Experiments Background*

Let *χ* be the compact design space and *y*(*t*) the response variable, assumed as a random variable of a parametric distribution with a vector of parameters *θ*. We assume that the observations *y*(*t*) verify the following model:

$$y(t) = \eta(t, \theta) + \varepsilon(t) \tag{4}$$

The model exhibits a random error *ε* and the function *η* is called the response surface and it is a partially known function, that is, it is within a parametric set of functions where the parameters *<sup>θ</sup><sup>T</sup>* = (*θ*1, ... , *<sup>θ</sup>k*) <sup>∈</sup> <sup>R</sup>*<sup>k</sup>* are unknown and their specification determines *<sup>η</sup>* totally. The response variable depends on the variable *t*, which we will call the independent or design variable, which can be freely chosen by the experimenter.

We will say that an exact design, *ξ<sup>N</sup>* is a collection of *N* different support points, *t*1, *t*2, ... , *tN*, where there could be repetitions. An approximate design, *ξ*, is a collection of *<sup>M</sup>* different support points, *<sup>t</sup>*1, *<sup>t</sup>*2, ... , *tM*, and weights, *<sup>ξ</sup>*(*t*1), *<sup>ξ</sup>*(*t*2), ... , *<sup>ξ</sup>*(*tM*), that defines a discrete probability at *<sup>χ</sup>*. If the design *<sup>ξ</sup>* has weight *<sup>ξ</sup>*(*ti*) at *ti*, *<sup>i</sup>* = 1, ... , *<sup>M</sup>* and the total number of observations made is *<sup>N</sup>*, approximately *<sup>N</sup>ξ*(*ti*) observations will be taken at *ti*, *i* = 1, ... , *M*. By choosing a sufficiently large *N*, it is possible to approximate a design with these characteristics to an exact design. The approximations will be better the larger the *N* [22].

One of the most important tools in the optimal experimental design is the FIM induced by *ξ*. The variance-covariance matrix *Var*[ˆ *θ*] of the maximum likelihood estimator ˆ *θ* can be asymptotically approximated by the inverse of the FIM, *M*−1(*ξ*), and under certain regularity conditions it reaches the lower limit of Cramer-Rao [23]. Once the estimator has been specified, we can consider the optimal design problem, as the selection of predictors that somehow lead to the minimization of the variance-covariance matrix *Var*[ˆ *<sup>θ</sup>*] <sup>≈</sup> *<sup>M</sup>*−1(*ξ*), or equivalently to the maximization of *M*(*ξ*). Therefore, the objective of the optimal design

of experiments will be to find the best of the designs that maximizes *M*(*ξ*) or minimizes *M*−1(*ξ*) in some sense.

Optimal design theory was initially developed for linear models. However, in practice it is common to find experiments where the response is explained from non-linear models in the parameters to be estimated, as this case study. The most common method for analysing data from a non-linear model is based on the use of the linear Taylor series approximation of the response surface.

$$\eta(t,\theta) \approx \eta(t,\theta^{(0)}) + [\nabla\_{\theta}\eta(\mathbf{x},\theta)|\_{\theta^{(0)}}]^T(\theta - \theta^{(0)})\tag{5}$$

where <sup>∇</sup>*<sup>θ</sup>* denotes the gradient with respect *<sup>θ</sup>* and being *<sup>θ</sup>*(0) a prior value of *<sup>θ</sup>*, for example an estimate calculated on the basis of an initial experiment or suitable values found in the literature. Thus, the FIM will depend on the values of the unknown parameters and, consequently, the optimal designs will also depend on these values,

$$M(\xi,\theta) = \int\_{\overline{\xi}} \nabla\_{\theta} \eta(t,\theta) \nabla\_{\theta}^{T} \eta(t,\theta) \xi(dt) \tag{6}$$

#### *2.2. Electrical Circuit Model*

Let us consider the EEC of Figure 1A. It comprises a resistance *R*<sup>1</sup> in series with the parallel combination of a capacitor *C* and a resistance *R*2. Currents and voltages have been labelled in each element. *i*(*t*) and *v*(*t*) are the current input and the voltage output, respectively. *iC*(*t*) and *iR*<sup>2</sup> (*t*) are the currents flowing through the capacitor *<sup>C</sup>* and the resistance *<sup>R</sup>*2, respectively. *vC*(*t*) is the voltage across *<sup>C</sup>* (the same than that of *<sup>R</sup>*2).

**Figure 1.** Electrical circuit model (**A**) and Nyquist plot (**B**).

The ODE of the EEC shown in Figure 1A is written as (refer to Appendix A)

$$\frac{d\upsilon(t)}{dt} + a\_0 \upsilon(t) = b\_1 \frac{d\dot{\imath}(t)}{dt} + b\_0 \dot{\imath}(t), \qquad t \in [0, T] \tag{7}$$

Let us define *<sup>θ</sup><sup>T</sup>* = (*b*1, *<sup>b</sup>*0, *<sup>a</sup>*0) for *<sup>b</sup>*1, *<sup>b</sup>*0, *<sup>a</sup>*<sup>0</sup> <sup>∈</sup> <sup>R</sup> and *<sup>k</sup>* <sup>=</sup> 3 the unknown vector parameters with

$$\begin{cases} b\_1 = R\_1 \\ b\_0 = \frac{R\_1 + R\_2}{R\_2 C} \\ a\_0 = \frac{1}{R\_2 C} \end{cases} \tag{8}$$

As we saw in Section 1, the impedance from Equation (7) is obtained as

$$Z(s) = \frac{b\_1s + b\_0}{s + a\_0} \tag{9}$$

For physical frequencies *s* = *jω*, Equation (9) is written as

$$Z(\omega) = \frac{b\_1 j \omega + b\_0}{j \omega + a\_0} \tag{10}$$

where *<sup>ω</sup>* <sup>∈</sup> (0, <sup>∞</sup>).

Figure 1B shows the Nyquist plot (−*X*(*ω*) vs. *<sup>R</sup>*(*ω*) in the complex plane). Equation (10) sketches a semicircle, which intersects the real axis at *<sup>R</sup>*<sup>1</sup> and *<sup>R</sup>*<sup>1</sup> + *<sup>R</sup>*2, at the frequencies *ω* = ∞ and *ω* = 0, respectively (*ω* is positive and increases counterclockwise). The maximum of <sup>−</sup>*X*(*ω*) vs. *<sup>R</sup>*(*ω*) is reached at *<sup>ω</sup>*max <sup>=</sup> <sup>1</sup> *<sup>R</sup>*2*<sup>C</sup>* .

Now, let the parameter *s* (refer to Equation (9)) be considered the differential operator, that is, *p* = *d*/*dt* (see above). The unknown parameters *θ* must be estimated from measures of the form:

$$v(t) = Z(p)i(t) + \varepsilon(t) \tag{11}$$

where *Z*(*p*) is a new differential operator [6] and *ε*(*t*) is a random error *N*(0, *σ*2) that includes both the errors of the tests and the specification of the model. The Fourier transform of Equation (11) is given by

$$V(\omega) = Z(\omega)I(\omega) + E(\omega) \tag{12}$$

where *<sup>V</sup>*(*ω*) = <sup>F</sup>[*v*](*ω*) <sup>∈</sup> <sup>C</sup>, *<sup>I</sup>*(*ω*) = <sup>F</sup>[*i*](*ω*) <sup>∈</sup> <sup>C</sup> and *<sup>E</sup>*(*ω*) = <sup>F</sup>[*ε*](*ω*) <sup>∈</sup> <sup>C</sup>, with <sup>F</sup> the Fourier transform.

It should be mentioned that the impedance of (9) can also be obtained by considering the impedance of each of *R*1, *R*<sup>2</sup> and *C*, that is, *R*1, *R*<sup>2</sup> and <sup>1</sup> *sC* , respectively, resulting in *<sup>Z</sup>*(*s*) = *<sup>R</sup>*<sup>1</sup> + (*R*<sup>2</sup> <sup>1</sup> *sC* ), denotes parallel, which gives Equation (9).

#### *2.3. The FIM for the Electrical Circuit Model*

In this section, we will consider that the input signal *i*(*t*) is a time stationary process and we will go from the time domain to the frequency domain in order to apply the optimal design theory for this type of LTI SISO (single-input, single-output) systems. In this domain, the properties for the approximation of the information matrices used in the theory of the optimal design of experiments, given by Kiefer-Wolfowitz [24], Karlin-Studden [25] and Fedorov [26], are preserved. According to Mehra [9], if a SISO system is LTI, the number of data points *N* is large and *i*(*t*) is a stationary process, the optimal inputs can be calculated much more efficiently using frequency domain techniques.

Let Equation (10) be the impedance model to be estimated, where *θ* is the vector of parameters of the impedance model defined in (8). Then,

$$\nabla\_{\theta}Z(\omega) = \begin{pmatrix} \frac{j\omega}{j\omega + a\_0} \\ \frac{1}{j\omega + a\_0} \\ -\frac{j(b\_1\omega - b\_0j)}{(j\omega + a\_0)^2} \end{pmatrix} \tag{13}$$

By Parseval's Theorem, we find that, for a continuous time model (10), the FIM in the frequency domain is given by (refer to Appendix B)

$$M(\Phi) = \frac{1}{2\pi} \int\_{-\infty}^{+\infty} \tilde{M}(\omega, \theta) \Phi(\omega) d\omega \tag{14}$$

where Φ is the power spectral density of *i*(*t*) and *M*˜ (*ω*, *θ*) is the single-frequency FIM, which defines the information associated with the model (10) with respect to *θ*, defined as

$$\tilde{M}(\omega,\theta) = \Re\left\{\nabla\_{\theta}Z(\omega) \cdot \nabla\_{\theta}Z(\omega)^{H}\right\} \tag{15}$$

where the superscript *H* denotes the conjugate transpose operator and {*z*} is the real part of the number *<sup>z</sup>* <sup>∈</sup> <sup>C</sup>. From the model (10), the single-frequency FIM as defined in (15) is

$$
\tilde{M}(\omega,\theta) = \begin{pmatrix}
\frac{\omega^2}{a\_0^2 + \omega^2} & 0 & \frac{\omega^2(b\_0 - b\_1 a\_0)}{(a\_0^2 + \omega^2)^2} \\
0 & \frac{1}{a\_0^2 + \omega^2} & -\frac{b\_1 \omega^2 + b\_0 a\_0}{(a\_0^2 + \omega^2)^2} \\
\frac{\omega^2(b\_0 - b\_1 a\_0)}{(a\_0^2 + \omega^2)^2} & -\frac{b\_1 \omega^2 + b\_0 a\_0}{(a\_0^2 + \omega^2)^2} & \frac{b\_1^2 \omega^2 + b\_0^2}{(a\_0^2 + \omega^2)^2}
\end{pmatrix} \tag{16}
$$

Note that in the frequency domain, the role of the probability measure is played by the power spectral density function Φ. Therefore, Φ is the design measure (or design), *ξ*, that we call in the context of design of experiments.

We define by Ξ the set of non-decreasing, continuous design measures of bounded variation, corresponding to input power spectral density functions, that is,

$$\Xi = \left\{ \Phi \; : \; \int\_{\omega\_I}^{\omega\_E} \Phi(\omega) d\omega = 1 \right\} \tag{17}$$

The design can be limited to a finite range [*ωI*, *<sup>ω</sup>E*], where *<sup>ω</sup><sup>I</sup>* and *<sup>ω</sup><sup>E</sup>* are the lowest and highest angular frequencies, respectively, at which impedance has been measured. This frequency range will be the design space *χ*.

#### **3. Construction of D-Optimal Input Signals**

Among a set of designs, it is not easy to decide which is "the best" of them. Therefore, we need to choose a criterion, a scalar measure of the FIM, <sup>Ψ</sup> : <sup>M</sup>(Ξ) −→ <sup>R</sup> ∪ {∞}, called optimality criterion, that helps us to find the best design. The choice will depend on the interests sought by conducting the experiment. A design that minimizes Ψ over all the designs on *χ* is called an optimum design, that is

$$\Phi^\* = \arg\min\_{\Phi \in \Sigma} \Psi[M(\Phi)] \tag{18}$$

The most popular criterion is D-optimality. This criterion consists of minimizing the volume of the confidence ellipsoid of estimators of the parameters of the model, i.e., it maximizes the determinant of the FIM. For a given *<sup>θ</sup>*(0), D-optimality is defined by the criterion function:

$$\Psi\_D[M(\Phi)] = [\det M(\Phi)]^{-\frac{1}{k}} \tag{19}$$

When an optimal design is sought among all approximate designs on the design space and the design criterion is convex, it can be checked the optimality of a particular design using the celebrated General Equivalence Theorem (GET) [24]. If the criterion function, Ψ, is differentiable the GET has a friendly version. Thus, a design Φ∗ *<sup>D</sup>* is Ψ*D*-optimum if and only if ∀*ω* ∈ *χ*, the dispersion function

$$d(\omega, \Phi\_D^\*) = \nabla\_\theta Z(\omega)^H M^{-1}(\Phi\_D^\*) \nabla\_\theta Z(\omega) = \text{Tr}\left[M^{-1}(\Phi\_D^\*) \tilde{M}(\omega, \theta)\right] \tag{20}$$

achieves its maximum values at the support points. In particular, for the D-optimum criterion this condition proves the optimality of Φ∗ *<sup>D</sup>* if and only if

$$d(\omega, \Phi\_D^\*) \le k, \ \forall \omega \in \chi \tag{21}$$

being *k* the number of the parameters of the model.

The GET allows us to check the optimality of a design, but it does not tell us how to find it. However, this theorem is useful for building efficient algorithms that allow computing optimal designs.

Another fundamental tool commonly used by the algorithmic techniques is the efficiency of a design. This value is interpretated as the goodness of a design and is defined as

*Eff*Ψ(Φ) = Ψ[*M*(Φ)]/Ψ[*M*(Φ∗)], being Φ<sup>∗</sup> the Ψ-optimal design. Thus the D-efficiency of a design Φ*<sup>D</sup>* is computed as:

$$Eff\_D(\Phi\_D) = \left[ \frac{\det M(\Phi\_D)}{\det M(\Phi\_D^\*)} \right]^{-\frac{1}{k}} \tag{22}$$

If the criteron function has a homogeneity property, as in this case, there is a practical statistical interpretation. Thus, if a design has 50 % efficiency, then half of the observations with the optimal design will produce the same results with respect to the criterion function. Although this quantity cannot be calculated when the optimal design is unknown, it is possible to obtain lower limits of efficiency to make a stopping rule. An important bound for the D-optimization criterion is the one proposed by Atwood [27]

$$Eff\_D(\Phi\_D) \ge \frac{k}{\max\_{\omega \in \chi} d(\omega, \Phi\_D)}\tag{23}$$

#### *Modified Algorithms*

The search for analytical solutions for the problem of construction of optimal designs turns out to be a difficult task and, in most of the real problems, it is not possible to calculate analytically the designs under a certain criterion. There is a rich literature on numerical computational algorithms proposed to obtain optimal designs under different scenarios. In exact design problems, we do not have a convex optimization problem in general and, so, finding optimal design is not an easy task because it is a discrete optimization problem and there is no general analytical tool for confirming whether an exact design is optimal or not. There are several numerical algorithms for finding optimal exact designs based on exchange methods. In this paper, we proposed the Atkinson and Donev KL exchange algorithm [28] to calculate exact D-optimal designs. The procedure starts with a non-singular random design. Then, two sets of points are constructed from the dispersion function. One with *K* "least promising" support points of the current design, for which exchanges are attempted. The other one with *L* "most promising" candidate design points. Finally, it adds and deletes observation points that lead to the greatest increase in the determinant of the FIM, under the standard constraint of the required number of points and the maximum execution time of the algorithm. Each exchange of improvement is executed immediately. Because this is only a heuristic, it cannot be guaranteed to reach the global optimum, so it is convenient to run the algorithm several times. The end result is the best design found in all runs.

For optimal approximate design problems, there are many analytic methods to construct them. Due to the performance and the flexibility to be applied in a broad range of problem structures and sizes, the randomized exchange algorithm (REX) [29] has been considered in this work. The REX method, which is a simple batch-randomized exchange algorithm, can be viewed as an efficient extension or combination of both the vertex exchange method for D-optimality of Böhning [30] and the KL exchange algorithm. The procedure begins with a random design of non-singular points and their respective proportions. Then, a batch consisting of a pair of points is repeatedly selected: a random support point of the current design, *ωk*, and a random design point, *ωl*, where optimal ratio exchanges between these pairs are made. The batch selection depends on the dispersion function. The key to making the optimal random exchanges of proportions between the selected batch is the optimal step length value *α*∗ *k*,*l* (Φ*D*) (see Appendix A in [29]), which gives the value to be added or subtracted from the proportions of each observation point. Finally, it adds and deletes observation points that lead to the largest increase in the efficiency bound. This algorithm has a standard threshold constraint on the minimum design efficiency to stop the computation. The limit used as the stopping rule is the Atwood limit (23).

These algorithms were adapted to involve complex variables. The matrix of Equation (13) is given as input for both algorithms. We develop a complex type null matrix and then it is completed by creating a loop, where each row of the complex null matrix is replaced by the components of (13). Each row corresponds to a value of *ω* within *χ* previously defined. Next, we adapt the definition of the single frequency FIM (see Equation (15)) and

the dispersion function (20). Modified codes to adapt the algorithms to this case study can be found in Appendix C.

The general procedure followed by these algorithms can be summarized as:

#### **4. Real Applications**

As mentioned earlier, biological tissues, cells, and ETIs can be characterized from their electrical behaviour. Figure 2 shows a portion of tissue (group of similar cells) with an implanted electrode. In general, the passive electrical behaviour of the tissue and the ETI involve electrical capacitance (capacity to store charges) and resistance (ability to oppose dc-current flow). Figure 2 also shows several basic EECs [1–4]. Intra- (pink colour) and extracellular (blue colour) media are resistive parts (electrolyte solutions), and thus are modelled by the resistances *RI* (Figure 2A) and *RE* (Figure 2B), respectively. The passive electrical behaviour of a cell membrane (green colour) is described by a parallel combination of a cell membrane capacitance *CM* and the membrane resistance *RM* (see Figure 2C). A basic model of ETI is shown in the EEC of Figure 2D, that is, a double layer capacitance *CDL* in parallel with the polarization resistance *RP* [1,3]. The EEC shown in Figure 1A involves those of Figure 2C,D with an additional resistance *R*1. Moreover, the impedance of the EEC of Figure 1A is equivalent to that proposed by Cole when the biological tissue involves an ideal capacitance [2,4].

It should be mentioned that the optimal design of experiments is of great interest for two main reasons: to obtain an optimal characterization of the process and to minimize the measurement acquisition time.

**Figure 2.** Biological tissue and basic EECs for intra- (**A**), extracellular media (**B**), cell membrane (**C**), and electrode-tissue interface (**D**).

#### *4.1. Methods*

The modified algorithms were applied to the circuit shown in Figure 1A. Let *<sup>R</sup>*<sup>1</sup> = <sup>100</sup> <sup>Ω</sup>, *<sup>R</sup>*<sup>2</sup> = 106 <sup>Ω</sup> and *<sup>C</sup>* = <sup>10</sup>−<sup>6</sup> <sup>F</sup> the nominal values of the parameters and *<sup>χ</sup>* = [2*π*10−1, 2*π*10] rad/s. For numerical study, we consider the frequency design space to be a set of grid points spread equidistantly by 0.01.

We have used the RStudio 1.3.1093 program to obtain the optimal designs. The Autolab PGSTAT204 potentiostat/galvanostat (Figure 3A), equipped with the FRA32M module, was used to perform the impedance measurements. The equipment, controlled by a computer and NOVA electrochemistry software, is connected to a dummy cell containing the circuit described above to perform the test (Figure 3B). We used sinusoidal current signals of 10−<sup>6</sup> A peak amplitude.

**Figure 3.** Autolab PGSTAT potentiostat/galvanostat equipment (**A**) connected to the dummy cell (**B**) containing the EEC.

The objective of this experimental test was to compare the efficiency of the classical design commonly used by experimenters with the D-optimal designs obtained in the adapted KL and REX algorithms. The classic design is based on measuring impedance at certain frequencies equally spaced on a logarithmic scale. In particular

$$\Phi^{\rm classical} = \{2\pi 10^{-1}, 2\pi 10^{-0.8}, 2\pi 10^{-0.6}, \dots, 2\pi 10^{0.8}, 2\pi 10\} \tag{24}$$

where *N* = 11. For comparison purposes, the adapted KL algorithm was executed considering the same number of support points as the classic design (*N* = 11). In the case of the REX algorithm, from a theoretical point of view, it is not necessary to set *N*. The algorithm provides points and the proportions to be measured at those points. In order to be able to compare the design obtained, the proportion was rounded to the nearest integer, that is, if the design puts weight *pi* on *<sup>ω</sup>i*, *<sup>i</sup>* = 1, ... , *<sup>M</sup>* and the total number of observations we want is *<sup>N</sup>* <sup>=</sup> 11, then approximately *<sup>N</sup>* · *pi* observations will be taken at *<sup>ω</sup>i*, *<sup>i</sup>* <sup>=</sup> 1, . . . , *<sup>M</sup>*.

#### *4.2. Results and Discussion*

Figure 4 shows the Nyquist plot for the EEC of Figure 3B. Experimental impedance data obtained for the frequencies of (24) draw a semicircle like that of Figure 1B. Subsequently, experimental data have been fitted to the EEC of Figure 1A and the values of its parameters are given in Table 1.

In this experimental test, there was a prior knowledge of the EEC parameters with a certain tolerance. In this case, the following initial values were considered to execute the algorithms:

$$b\_1 = 10, \quad b\_0 = 50500, \quad a\_0 = 50 \tag{25}$$

where *<sup>R</sup>*<sup>1</sup> <sup>=</sup> <sup>10</sup> <sup>Ω</sup>, *<sup>R</sup>*<sup>2</sup> <sup>=</sup> <sup>1000</sup> <sup>Ω</sup> and *<sup>C</sup>* <sup>=</sup> <sup>20</sup> <sup>×</sup> <sup>10</sup>−<sup>6</sup> F. To avoid dependence on the value of these parameters, a sequential approach was applied with the estimates obtained in the previous step. The stabilization of the process was fast (Table 1). Finally, the optimal design obtained with the two algorithms was the following:

$$
\Phi\_D^\* = \left\{ \begin{array}{cc} 2\pi 0.1 & 2\pi 10 \\ 7 & 4 \end{array} \right\} \tag{26}
$$

**Figure 4.** Nyquist plot obtained from the classic design.

**Table 1.** EEC parameters values estimated using the different optimal designs.

The dispersion function (20), shown in Figure 5, reaches the maximums at the support points (that is, the extremes of the interval) and remains less than 3 (the number of parameters) throughout the design space, so the design obtained is D-optimal.

**Figure 5.** Dispersion function for the D-optimal design for the EEC shown in Figure 1A. The abscissa axis represents the logarithm of the frequency *f* (*f* = *ω*/2*π*).

Figure 6 illustrates the support points of the D-optimal designs obtained with *N* = 11 after executing the algorithms and the classic design. Table 2 shows the values of the D-optimal criterion (19) for the designs obtained from REX, KL and classical. It is observed that the values obtained with the algorithms are greater than the value obtained with the classical design, therefore, the value of the D-optimal criterion has been maximized. Optimal designs are an interesting tool to measure the value of an experimental design, through efficiency. This efficiency is the percentage of observations that the optimal design would need to achieve the precision of the experimental design being compared. The efficiency (22) of the classic design is also presented in Table 2. This shows how the same information can be obtained by reducing the number of experimental conditions by using optimal inputs.

**Figure 6.** Support points of the D-optimal designs obtained by the adapted algorithms and the classic design. The values on the points obtained with the algorithms (7 and 4) are the number of replicas of each point.

**Table 2.** Criterion values.

Finally, a sensitivity study has been performed by analysing the fluctuations of the efficiency values against the initial variations of the parameters. In particular, each initial parameter was modified by ±20 % of the real values. The efficiencies of the designs obtained have been calculated in relation to the 3<sup>3</sup> = 27 possible designs. For all cases, as can be seen from Figure 7, the efficiency for both REX and the KL was never lower than 99.21 % and higher than 99.99 % for most variations.

**Figure 7.** Efficiencies of the 27 possible designs.

#### **5. Conclusions**

This paper expands the existing tools to create optimal designs regarding SISO LTI systems in the context of electrical impedance measurements, based on the frequency domain description of an input signal. We have considered the most used optimization criterion, D-optimality. Two algorithms (KL and REX) were modified to calculate Doptimal designs and estimate the parameters of the EEC describing the electrical behaviour of bioelectrodes, cell membranes, or biological tissues. We have checked that the efficiency remains high for different initial values of the parameters. The estimation of the parameters has also been discussed. The application of a sequential approximation with the estimates obtained by means of least squares stabilize the process quickly. We have calculated the D-optimal designs with the same number of points as those of the classical design to obtain an efficiency comparison. It has been seen that optimal designs computed may save 35% of the observations to get the same results as the classic design.

The generation of new results for the design of optimal input signals in impedance models is very useful both in the characterization of biological tissues and the electrode– tissue interface. Specifically, the KL exchange and random exchange REX algorithms that we have adapted are a good option to be used by experimenters to obtain a good estimate of the parameters with lower time-consuming.

**Author Contributions:** Conceptualization, R.M.M. and J.-L.P.S.; methodology, À.S.B., R.M.M. and J.-L.P.S.; software, À.S.B. and R.M.M.; writing—original draft preparation, À.S.B. and R.M.M.; writing—review and editing, À.S.B., J.-L.P.S. and R.M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Ministerio de Economía, Industria y Competitividad-Agencia Estatal de Investigación (Project MTM2016-80539-C2-1-R), the Consejería de Educación, Cultura y Deportes of Junta de Comunidades de Castilla-La Mancha (Project SBPLY/17/180501/000380), and the European Regional Development Fund (FONDOS FEDER).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Obtaining the Differential Equation**

Kirchhoff's current law provides

$$i(t) = i\_{\mathbb{C}}(t) + i\_{\mathbb{R}\_2}(t) \tag{A1}$$

Applying Kirchhoff's voltage law, we obtain

$$
\upsilon\_{\mathbb{C}}(t) = \upsilon(t) - \upsilon\_{R\_1}(t) \tag{A2}
$$

As shown in Figure 1A, *C* and *R*<sup>2</sup> are in parallel, and therefore have the same voltage *vC*(*t*). Using the characteristic equations of the elements

$$
v\_{\mathbb{C}}(t) = \mathcal{R}\_2 i\_{\mathbb{R}\_2}(t) \tag{A3}$$

$$
\sigma\_{R\_1}(t) = R\_1 i(t) \tag{A4}
$$

$$i\_C(t) = \mathbb{C} \frac{dv\_C(t)}{dt} \tag{A5}$$

Using the Equations (A3) and (A5) in Equation (A1), we obtain

$$i(t) = \mathcal{C}\frac{dv\_{\mathcal{C}}(t)}{dt} + \frac{v\_{\mathcal{C}}(t)}{R\_2} \tag{A6}$$

And, using Equations (A2) and (A4) in Equation (A6), we obtain

$$\dot{v}(t) = \mathcal{C}\frac{dv(t)}{dt} - \mathcal{C}R\_1\frac{di(t)}{dt} + \frac{v(t)}{R\_2} - \frac{R\_1}{R\_2}i(t) \tag{A7}$$

Finally, dividing both members by *C* and rearranging we arrive at the ODE:

$$\frac{d\upsilon(t)}{dt} + \frac{1}{R\_2 \mathcal{C}} \upsilon(t) = R\_1 \frac{d\dot{\imath}(t)}{dt} + \frac{R\_1 + R\_2}{R\_2 \mathcal{C}} \dot{\imath}(t), \quad t \in [0, T] \tag{A8}$$

#### **Appendix B. Fisher Information Matrix of a LTI System**

*Appendix B.1. Time Domain*

The impulse response fully characterizes LTI systems and allows the output of an LTI system to be calculated against any input. As indicated above, an LTI circuit, with input *i*(*t*) and output *v*(*t*), can be described by an ODE with constant coefficients (see Equation (1)). The output *v*(*t*) of the continuous LTI system is given by a convolution integral:

$$v(t) = z(t) \* i(t) = \int\_{-\infty}^{+\infty} z(t - \tau)i(\tau)d\tau \tag{A9}$$

where *z*(*t*) is the voltage response to a unit-impulse current at the input. Therefore, taking into account that the FIM for a non-linear model has the form

$$\begin{split} M(i(t),\theta) &= \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \nabla\_{\theta} z(t-\tau,\theta) \nabla\_{\theta}^{T} z(t-\nu,\theta) i(\tau) i(\nu) d\tau \, d\nu \, dt \\ &= \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \nabla\_{\theta} z(\tau,\theta) \nabla\_{\theta}^{T} z(\nu,\theta) \underbrace{i(t-\tau) i(t-\nu)}\_{(\star)} d\tau \, d\nu \, dt \end{split} \tag{A10}$$

The function () will be needed for the next subsection.

#### *Appendix B.2. Frequency Domain*

Assuming that *i*(*t*) is a time stationary process (the value of the covariance between two periods depends only on the distance between these two periods of time) and applying the limit to the function (), we obtain the next

$$\lim\_{T \to \infty} \frac{1}{2T} \int\_{-T}^{T} i(t - \tau)i(t - \nu)dt = \lim\_{T \to \infty} \frac{1}{2T} \int\_{-T}^{T} i(t)i(t + \tau - \nu) = R\_i(\tau - \nu) \tag{A11}$$

where the time period is *<sup>t</sup>* <sup>−</sup> *<sup>ν</sup>* <sup>−</sup> (*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*) = *<sup>τ</sup>* <sup>−</sup> *<sup>ν</sup>* and *Ri* is the autocorrelation function of the signal *i*(*t*).

The normalized FIM is defined as

$$\begin{split} M(R\_i(\tau)) &= \lim\_{T \to \infty} \frac{1}{2T} M(i(t), \theta) \\ &= \lim\_{T \to \infty} \frac{1}{2T} \int\_{-T}^{T} \int\_{-T}^{T} \int\_{-T}^{T} \nabla\_{\theta} z(\tau, \theta) \nabla\_{\theta}^{\top} z(\nu, \theta) i(t - \tau) i(t - \nu) d\tau \, d\nu \, dt \\ &= \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \nabla\_{\theta} z(\tau, \theta) \nabla\_{\theta}^{\top} z(\nu, \theta) \lim\_{T \to \infty} \frac{1}{2T} \left[ \int\_{-T}^{T} i(t - \tau) i(t - \nu) dt \right] d\tau d\nu \\ &= \int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \nabla\_{\theta} z(\tau, \theta) \nabla\_{\theta}^{\top} z(\nu, \theta) R\_i(\tau - \nu) \, d\tau \, d\nu \end{split} \tag{A12}$$

The Wiener–Khinchin theorem establishes a relationship between the autocorrelation function of a signal and its power spectral density function, that is, the Fourier transform of the autocorrelation function *Ri*(*τ*) of a time stationary process is the function spectral density <sup>Φ</sup>(*ω*) (F[*Ri*](*ω*) = <sup>Φ</sup>(*ω*)). Therefore, by Parseval's Theorem, we find that the FIM in the frequency domain is given by

$$\begin{split} M(\mathcal{F}[R\_{i}]) \underset{\text{Notation}}{\equiv} M(\Phi) &= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \mathcal{F} \Big[\nabla\_{\theta} z(\tau, \theta)\Big] \mathcal{F} \Big[\nabla\_{\theta} z(\nu, \theta)\Big] \mathcal{F} (R\_{i}(\tau - \nu)) d\omega \\ &= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \Re\Big\{\nabla\_{\theta} \mathcal{F} \Big[z(\tau, \theta)\Big] \nabla\_{\theta} \mathcal{F} \Big[z(\nu, \theta)\Big] \right\} \Phi(\omega) d\omega \\ &= \frac{1}{2\pi} \int\_{-\infty}^{\infty} \Re\Big\{\nabla\_{\theta} \mathcal{Z}(\omega) \nabla\_{\theta} \mathcal{Z}(\omega)^{H}\Big\} \Phi(\omega) d\omega \end{split} \tag{A13}$$

where {*z*} is the real part of the number *<sup>z</sup>* <sup>∈</sup> <sup>C</sup>, <sup>∇</sup>*<sup>θ</sup>* <sup>=</sup> *<sup>∂</sup> ∂θ* and the superscript *H* denotes the conjugate transpose operator. Indeed, the Fourier integral of *z*(*τ*, *θ*)*z*(*ν*, *θ*) is the function *Z*(*ω*)*Z*(*ω*)*H*.

#### **Appendix C. Modified Algorithm Codes**

```
space <−seq ( 0 . 6 2 , 6 2 . 8 3 , by =0.01)
n <−len g th ( space )
m <−3
MATRIX_F <− ma trix ( 0 ^ 1 i , n , m)
b1 <−100
b0 <−1000100
a0 <−1
j <−1
for ( i in space )
{
f <− f u n c ti o n ( x )
c ( complex ( re al =0 ,imaginary=x )/complex ( re al=a0 , imaginary=x ) ,
1/complex ( re al=a0 , imaginary=x ) ,
complex ( re al=b0 , imaginary=b1 * x)/
( complex ( r e a l =x , imaginary=−a0 ) ) ^ 2 )
```
M\_row <− f ( i ) MATRIX\_F [ j , ] <− M\_row j <− j +1 } *Appendix C.1. REX* Modification of the FIM: M <−( t ( s q r t (w. supp )\* Conj ( Fx . supp ) ) %\*% ( s q r t (w. supp )\* Fx . supp ) ) M <−Re (M) Modification of the dispersion function: Mp <− ma trix (NA, n , 1 ) for ( i in 1:n) { Mpp <− Fx [ i , ] %\*% t ( Conj ( Fx [ i , ] ) ) Mpp <− Re (Mpp) Z <− s ol ve (M+E ) %\*% Mpp t r a z a <−sum ( diag (Z ) ) Mp[ i , ] <− t r a z a } *Appendix C.2. KL* Modification of the FIM: M1 <− t ( (w %\*% one ) \* F) %\*% Conj (F ) M2 <− Re (M1) Modification of the dispersion function: H <− Conj ( F ) %\*% s ol ve (M2 + E )

H2 <− H \* F d . fun <− apply ( Re (H2 ) , 1 , sum )

#### **References**

### *Article* **Discriminative Convolutional Sparse Coding of ECG Signals for Automated Recognition of Cardiac Arrhythmias**

**Bing Zhang and Jizhong Liu \***

Nanchang Key Laboratory of Medical and Technology Research, School of Advanced Manufacturing, Nanchang University, Nanchang 330031, China

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Electrocardiogram (ECG) is a common and powerful tool for studying heart function and diagnosing several abnormal arrhythmias. In this paper, we present a novel classification model that combines the discriminative convolutional sparse coding (DCSC) framework with the linear support vector machine (LSVM) classification strategy. In the training phase, most existing convolutional sparse coding frameworks are unsupervised in the sense that label information is ignored in the convolutional filter training stage. In this work, we explicitly incorporate a label consistency constraint called "discriminative sparse-code error" into the objective function to learn discriminative dictionary filters for sparse coding. The learned dictionary filters encourage signals from the same class to have similar sparse codes, and signals from different classes to have dissimilar sparse codes. To reduce the computational complexity, we propose to perform a max-pooling operation on the sparse coefficients. Using LSVM as a classifier, we examine the performance of the proposed classification system on the MIT-BIH arrhythmia database in accordance with the AAMI EC57 standard. The experimental results show that the proposed DCSC + LSVM algorithm can obtain 99.32% classification accuracy for cardiac arrhythmia recognition.

**Keywords:** electrocardiogram signal; discriminative convolutional sparse coding; dictionary filter learning; linear SVM

**MSC:** 68U01

#### **1. Introduction**

Electrocardiogram (ECG) is used to record cardiac activity and detect different abnormalities in cardiac function and is a commonly used non-invasive tool for non-invasive diagnosis of cardiac arrhythmias. However, visual analysis is extremely limited and imprecise due to the large amount of information contained in the ECG, which may lead to misdiagnosis or inaccurate detection of arrhythmias. Therefore, computer-aided analysis helps doctors to detect cardiac arrhythmias quickly and efficiently.

The automatic arrhythmia detection system mainly includes feature extraction, feature selection, and classifier construction. ECG signal feature extraction techniques can be divided into time-based methods [1–3], frequency methods [4–6], and time-frequency techniques [7,8]. The time domain features mainly include heartbeat interval, duration parameters, and amplitude parameters. Due to the subtle changes in ECG amplitude and duration, time-based methods do not provide good discrimination [9,10]. Therefore, frequency methods using such as Fourier transform and power spectral density (PSD) and time-frequency methods using wavelet transform are proposed. However, the frequency method does not provide time information from the ECG signal. Time-frequency technology based on wavelet transform is widely used in time-frequency feature extraction of ECG signals. Before the time-frequency feature vectors extracted by wavelet transform are applied to the classifier, it is important to choose the best dimensionality reduction method. Dimensionality reduction methods for linear and non-linear transforms based on wavelet

**Citation:** Zhang, B.; Liu, J. Discriminative Convolutional Sparse Coding of ECG Signals for Automated Recognition of Cardiac Arrhythmias. *Mathematics* **2022**, *10*, 2874. https://doi.org/10.3390/ math10162874

Academic Editors: Carmen Lacave and Ana Isabel Molina

Received: 13 July 2022 Accepted: 9 August 2022 Published: 11 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

transforms have been proposed [7,11–13]. Martis et al. [7] transformed the ECG heartbeat using DWT and then applied independent component analysis (ICA) methods to extract features. The experimental results show that the features extracted using the ICA method combined with a probabilistic neural network (PNN) have better classification results. In the classification step, the most commonly used techniques in ECG signal classification are support vector machines (SVM) [1,14] and artificial neural networks [15].

In this paper, unlike most existing methods, we propose a novel classification model that combines the discriminative convolutional sparse coding framework with the linear support vector machines (LSVM) classification strategy.

The convolutional sparse coding (CSC) model assumes that the signal can be represented as a superposition of a few local filters, convolved with sparse feature maps. The CSC model handles the signal globally, and yet pursuit and dictionary learning are feasible due to the specific structure of the dictionary involved. CSC has been utilized for a variety of computer vision and pattern recognition tasks, such as inpainting [16], image separation [17], image fusion [18], object recognition [19], pedestrian detection [20], and tissue classification [21].

Although the CSC model has achieved state-of-the-art performance in many fields, most of the existing CSC frameworks are unsupervised and ignore label information in the convolutional filter training stage. Chen et al. [22] proposed a novel convolutional sparse coding classification (CSCC) approach, which introduces label information during the training process. The results show that dictionaries trained by convolution with label information can obtain more representative image information and achieve better classification performance. However, in the CSCC model, a set of dictionary filters are trained for each class of data, and then the dictionary filters from all classes build the final discriminative filter bank. This is not an efficient approach. At the same time, the intra-class and inter-class information are not considered, so the coding coefficients have small within-class scatter but large between-class scatter.

How to improve the classification performance of convolutional sparse coding has become a research direction in this field. One of the solutions is to obtain a discriminative dictionary filter bank through training so that the samples of different classes encoded by the dictionary filter bank are discriminative. Methods for learning a discriminative dictionary for sparse coding from training data have been recently proposed. Zhang and Li [23] proposed a discriminative KSVD dictionary learning algorithm. The algorithm incorporated the "classification error" term into the objective function, obtained the dictionary through the KSVD algorithm, and finally realized the classification task through a linear classifier. Yang et al. [24] learned a structured dictionary with class labels by adding Fisher's discriminant criterion during encoding. The method considers the intra-class and inter-class information in the coding process so that the coding coefficients have small within-class scatter but large between-class scatter. Following the work in [23], Jiang et al. [25] proposed to incorporate a "discriminative sparse code error" term into the objective function to enhance the discriminative power. It forces the signals from the same class to have very similar sparse representations, which results in good classification performance even using a simple linear classifier.

Obtaining a discriminative dictionary for classification by sparse coding has been successfully applied and achieved excellent performance. However, the discriminative dictionary obtained by sparse coding also has shortcomings, that is, it cannot capture the shifted local features in the sample, and at the same time, for high-dimensional signals, there is a curse of dimensionality. The dictionary filter bank obtained by convolutional sparse coding has shift invariance, local features at the sample translation position are extracted by convolution, and there is no curse of dimensionality.

However, how to learn discriminative convolutional dictionary filters with convolutional sparse representation and classification functions by supervised training is still an issue worth investigating. In this work, we propose a discriminative convolutional sparse coding (DCSC) model to learn discriminative dictionary filters for sparse coding of ECG

signals. Different from the CSCC algorithm, we explicitly incorporate a label consistency constraint called "discriminative sparse-code error" into the objective function, transform the DCSC model to the Fourier domain to reduce the computational cost [26], and optimize the DCSC model using the alternating direction method of the multiplier framework [27] (ADMM) algorithm. By adding label information to the training phase of the convolutional dictionary filters, the learned dictionary filters encourage signals from the same class to have similar sparse codes, and signals from different classes to have dissimilar sparse codes. To reduce the computational complexity, we use the max-pooling operation on the sparse coefficients. These pooled coefficients are then used as features and fed to an LSVM classifier for the ECG classification task.

Our contributions are as follows:

The rest of the paper is organized as follows: Section 2 presents the background of CSC, CDL, and LC-KSVD. Section 3 describes the proposed DCSC model. Section 4 describes the experimental results. The conclusion is given in Section 5.

#### **2. Literature Survey**

#### *2.1. Convolutional Sparse Coding*

The convolutional sparse coding (CSC) model assumes that a signal **<sup>s</sup>** ∈ *<sup>N</sup>* can be represented by the sum of *M* convolutions. These are built by feature maps {**x***m*}, each of length *N*, convolved with *M* small support filters {**d***m*} of length *n N*, ∗ denotes convolution, and the CSC model can be formulated as:

$$\underset{\{\mathbf{x}\_{m}\}}{\operatorname{argmin}} \frac{1}{2} \left\| \sum\_{m} \mathbf{d}\_{m} \ast \mathbf{x}\_{m} - \mathbf{s} \right\|\_{2}^{2} + \lambda \sum\_{m} \left\| \mathbf{x}\_{m} \right\|\_{1} \tag{1}$$

where *λ* > 0 is a regularization parameter. Given the filters, the above problem becomes the CSC pursuit task of finding the representations {**x***m*}.

#### *2.2. Convolutional Dictionary Learning*

A common approach for convolutional dictionary learning (CDL) [26] entails using a batch of *K* training data to optimize the filters and sparse coefficient maps. This problem can be formulated as follows:

$$\underset{\{\mathbf{d}\_{m}\}\_{1},\{\mathbf{x}\_{m,k}\}}{\operatorname{argmin}} \frac{1}{2} \sum\_{k} \left\| \sum\_{m} \mathbf{d}\_{m} \ast \mathbf{x}\_{m,k} - \mathbf{s}\_{k} \right\|\_{2}^{2} + \lambda \sum\_{k} \sum\_{m} \left\| \mathbf{x}\_{m,k} \right\|\_{1} \text{ such that } \left\| \mathbf{d}\_{m} \right\|\_{2} = 1 \,\forall m \tag{2}$$

where the constraint on the norms of dictionary filters {**d***m*} is required to avoid the scaling ambiguity between filters and coefficients. We denote the number of filters and the amount of training data by *M* and *K*, respectively.

The CDL problem is usually addressed by alternating optimization with respect to **x***m*,*<sup>k</sup>* and {**d***m*}. Several works have shown that solving (2) with respect to {**d***m*} can also be done effectively and efficiently using ADMM in the frequency domain [26,28].

#### *2.3. Label Consistent KSVD*

Jiang et al. [25] proposed a label consistent KSVD (LC-KSVD1) algorithm to learn a discriminative dictionary. The objective function is defined as follows:

$$<\mathbf{D}, \mathbf{A}, \mathbf{X} > = \arg\min\_{\mathbf{D}, \mathbf{A}, \mathbf{X}} ||\mathbf{S} - \mathbf{D}\mathbf{X}||\_2^2 + a ||\mathbf{Q} - \mathbf{A}\mathbf{X}||\_2^2 \quad \text{s.t.} \ ||\mathbf{x}\_i||\_0 < T \tag{3}$$

where ||**<sup>S</sup>** <sup>−</sup> **DX**||<sup>2</sup> <sup>2</sup> denotes the reconstruction error and **<sup>D</sup>** is the learned dictionary. ||**<sup>Q</sup>** <sup>−</sup> **AX**||<sup>2</sup> 2 represents the discriminative sparse-code error, **Q** are the "discriminative" sparse codes of input signals **S** used for classification, and **A** is a linear transformation matrix. *α* is the scalar that controls the relative contribution of the corresponding term. **X** are the sparse codes of input signals **S** and **x***<sup>i</sup>* is the column vector of **X**. *T* is a sparsity constraint factor (each signal has fewer than *T* non-zero items in its decomposition). By adding a discriminative sparse-code error constraint, the LC-KSVD1 model has good classification performance even with a simple linear classifier.

#### **3. The Proposed ECG Signal Classification System**

The proposed DCSC-based automatic arrhythmia identification system framework is shown in Figure 1. The proposed method consists of four stages, namely (1) in the discriminative convolutional sparse dictionary filters learning stage, we trained the DCSC model to obtain a discriminative dictionary filter. (2) In the sparse coding stage, we used the CSC model to obtain the discriminative sparse coefficients of the training and the test signals. (3) Pooling was performed on those sparse coefficients to reduce the large amount of data to an appropriate level. More importantly, it is used to obtain a compact representation of features that are invariant to local transformations. (4) In the ECG heartbeat testing stage, we used the pooled coefficients as features and fed them to an LSVM classifier for ECG classification.

**Figure 1.** Block diagram of the proposed ECG classification system.

#### *3.1. Discriminative Convolutional Sparse Dictionary Learning Model*

To improve the discriminative properties of the convolutional sparse coefficients, we incorporate the discriminative terms into the objective function during training. It can be formulated as the following optimization problem:

argmin {**d***m*},{**a***m*},{**x***k*,*m*} 1 <sup>2</sup> ∑ *k* ( ( (∑ *m* **d***<sup>m</sup>* ∗ **x***m*,*<sup>k</sup>* − **s***<sup>k</sup>* ( ( ( 2 2 + 1 <sup>2</sup> *α* ∑ *k* ( ( (∑ *m* **a***<sup>m</sup>* ∗ **x***m*,*<sup>k</sup>* − **q***<sup>k</sup>* ( ( ( 2 2 + *λ* ∑ *k* ∑ *m* **x***m*,*k*<sup>1</sup> such that **d***m*<sup>2</sup> <sup>=</sup> 1, **a***m*<sup>2</sup> <sup>=</sup> <sup>1</sup>∀*<sup>m</sup>* (4)

where **<sup>s</sup>***<sup>k</sup>* is the *<sup>k</sup>th* training data set, **<sup>s</sup>***<sup>k</sup>* ∈ *N*. *<sup>λ</sup>* <sup>&</sup>gt; 0 is a regularization parameter and *<sup>α</sup>* is the scalar that controls the relative contribution of the corresponding terms. {**d***m*} and {**a***m*} are a set of *<sup>M</sup>* dictionary filters, and **x***m*,*<sup>k</sup>* is a set of coefficient maps corresponding to the *mth* dictionary filter and the *kth* training data, each having the same size as **s***k*. **<sup>q</sup>***<sup>k</sup>* <sup>=</sup> [0, . . . , 1, 1, . . . , 0] *<sup>t</sup>* ∈ *<sup>N</sup>* is a "discriminative" sparse code corresponding to an input signal **s***k*, and **q***<sup>k</sup>* is the column vector of **Q**, whose size is the same as **s***k*. For example, assuming that **<sup>S</sup>** = [**s**0, **<sup>s</sup>**1, **<sup>s</sup>**2, **<sup>s</sup>**3], where **<sup>s</sup>**0, **<sup>s</sup>**<sup>1</sup> are from class 1, **<sup>s</sup>**2, **<sup>s</sup>**<sup>3</sup> are from class 2, and the length of the training data set *N* = 6, then **Q** can be defined as:

$$\mathbf{Q} = \begin{bmatrix} 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 1 & 1 \end{bmatrix} \tag{5}$$

The codes **x***m*,*<sup>k</sup>* and the dictionary filters {**d***m*} and {**a***m*} can then be efficiently updated in an alternating manner, as follows.

#### A. Convolutional sparse coding (CSC) step

Given the dictionary filters {**d***m*} and {**a***m*}, the codes can be updated by using ADMM in the frequency domain.

$$\underset{\{\mathbf{x}\_{m,k}\}}{\operatorname{argmin}} \frac{1}{2} \sum\_{k} \left\| \sum\_{m} \mathbf{d}\_{m} \ast \mathbf{x}\_{m,k} - \mathbf{s}\_{k} \right\|\_{2}^{2} + \frac{1}{2} a \sum\_{k} \left\| \sum\_{m} \mathbf{a}\_{m} \ast \mathbf{x}\_{m,k} - \mathbf{q}\_{k} \right\|\_{2}^{2} + \lambda \sum\_{k} \sum\_{m} \left\| \mathbf{x}\_{m,k} \right\|\_{1} \tag{6}$$

If we define the linear operators **<sup>D</sup>***<sup>m</sup>* and **<sup>A</sup>***<sup>m</sup>* such that **<sup>D</sup>***m***x***m*,*<sup>k</sup>* <sup>=</sup> **<sup>d</sup>***<sup>m</sup>* <sup>∗</sup> **<sup>x</sup>***m*,*k*, **<sup>A</sup>***m***x***m*,*<sup>k</sup>* <sup>=</sup> **a***<sup>m</sup>* ∗ **x***m*,*k*, and denote

**<sup>D</sup>** = (**D**<sup>0</sup> **<sup>D</sup>**<sup>1</sup> ...), **<sup>A</sup>** = (**A**<sup>0</sup> **<sup>A</sup>**<sup>1</sup> ...), **<sup>S</sup>** = (**s**<sup>0</sup> **<sup>s</sup>**<sup>1</sup> ...), **<sup>Q</sup>** = (**q**<sup>0</sup> **<sup>q</sup>**<sup>1</sup> ...), **<sup>X</sup>** <sup>=</sup> ⎛ ⎜⎝ **x**0,0 **x**0,0 ··· **x**1,0 **x**0,0 ··· . . . . . . ... ⎞ ⎟⎠ (7)

then we can rewrite Equation (6) as

$$\underset{\mathbf{X}}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{D}\mathbf{X} - \mathbf{S}\|\_{2}^{2} + \frac{1}{2} \kappa \|\mathbf{A}\mathbf{X} - \mathbf{Q}\|\_{2}^{2} + \lambda \|\mathbf{X}\|\_{1} \tag{8}$$

The ADMM algorithm and shrinkage/soft thresholding algorithm can be employed to solve Equation (8), as described in Appendix A.

B. Convolutional dictionary update (CDU) step

In developing the dictionary filter update, it is convenient to switch the index of the coefficient map from **x***m*,*<sup>k</sup>* to **x***k*,*<sup>m</sup>* . With the codes fixed, the dictionary filters {**d***m*} and {**a***m*} can be updated by solving the following optimization problems,

$$\begin{aligned} \underset{\{\mathbf{d}\_{m}\},\{\mathbf{a}\_{w}\},\{\mathbf{a}\_{k,w}\}}{\arg\min} \quad & \frac{1}{2} \sum\_{k} \left\| \sum\_{m} \mathbf{x}\_{k,m} \ast \mathbf{d}\_{m} - \mathbf{s}\_{k} \right\|\_{2}^{2} + \frac{1}{2} a \sum\_{k} \left\| \sum\_{m} \mathbf{x}\_{k,m} \ast \mathbf{a}\_{m} - \mathbf{q}\_{k} \right\|\_{2}^{2} \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \text{such that } \left\| \mathbf{d}\_{m} \right\|\_{2} = 1, \left\| \mathbf{a}\_{m} \right\|\_{2} = 1 \forall m \end{aligned} \tag{9}$$

Updating {**d***m*}

With the **x***k*,*<sup>m</sup>* and {**a***m*} fixed, the dictionary filters {**d***m*} can be updated by solving the following optimization problem,

$$\underset{\{\mathbf{d}\_m\}}{\operatorname{argmin}} \frac{1}{2} \sum\_{k} \left\| \sum\_{m} \mathbf{x}\_{k,m} \ast \mathbf{d}\_m - \mathbf{s}\_k \right\|\_2^2 \text{ such that } \|\mathbf{d}\_m\|\_2 = 1 \,\forall m \tag{10}$$