Experimental designs for a Benign Paroxysmal Positional Vertigo model

Background The pathology of the Benign Paroxysmal Positional Vertigo (BPPV) is detected by a clinician through maneuvers consisting of a series of consecutive head turns that trigger the symptoms of vertigo in patient. A statistical model based on a new maneuver has been developed in order to calculate the volume of endolymph displaced after the maneuver. Methods A simplification of the Navier‐Stokes problem from the fluids theory has been used to construct the model. In addition, the same cubic splines that are commonly used in kinematic control of robots were used to obtain an appropriate description of the different maneuvers. Then experimental designs were computed to obtain an optimal estimate of the model. Results D‐optimal and c‐optimal designs of experiments have been calculated. These experiments consist of a series of specific head turns of duration Δt and angle α that should be performed by the clinician on the patient. The experimental designs obtained indicate the duration and angle of the maneuver to be performed as well as the corresponding proportion of replicates. Thus, in the D‐optimal design for 100 experiments, the maneuver consisting of a positive 30° pitch from the upright position, followed by a positive 30° roll, both with a duration of one and a half seconds is repeated 47 times. Then the maneuver with 60° /6° pitch/roll during half a second is repeated 16 times and the maneuver 90° /90° pitch/roll during half a second is repeated 37 times. Other designs with significant differences are computed and compared. Conclusions A biomechanical model was derived to provide a quantitative basis for the detection of BPPV. The robustness study for the D‐optimal design, with respect to the choice of the nominal values of the parameters, shows high efficiencies for small variations and provides a guide to the researcher. Furthermore, c‐optimal designs give valuable assistance to check how efficient the D‐optimal design is for the estimation of each of the parameters. The experimental designs provided in this paper allow the physician to validate the model. The authors of the paper have held consultations with an ENT consultant in order to align the outline more closely to practical scenarios.


Introduction
First described by Bárány [1], Benign Paroxysmal Positional Vertigo (BPPV) is the most common vestibular disorder leading to vertigo. These vestibular symptoms are precipitated when the orientation of the head or body is changed relative to gravity, provoking http://www.tbiomed.com/content/10/1/21 brief periods (2-3 minutes) of vertigo, imbalance, and nausea. These changes can occur during daily activities such as lying down in bed or reaching up to retrieve an object from a high shelf. Benign Paroxysmal Positional Vertigo is commonly called top-shelf vertigo [2].
Contrary to what is widely believed, this type of vertigo is caused by a disorder in an organ of the inner ear, called the semicircular canal, which regulates balance. Figure 1 illustrates the three canals and how they are arranged in a similar way to the three cartesian axes. Each canal is filled with a fluid called endolymph and contains motion sensors within the fluids. At the base of each canal, the bony region of the canal has a dilated sac at one end called the ampulla. Within the ampulla is a mound of hair cells and supporting cells called crista ampullaris. These hair cells are composed of many cilia and are embedded in a gelatinous structure called cupula. As the head rotates, the duct moves, but the endolymph lags behind. This deflects the cupula and bends the cilia within. The bending of these cilia alters an electric signal that is transmitted to the brain, which sends this information to the eyes, provoking the corresponding vestibular movement which helps us keep our balance.
When does BPPV occur? Semicircular canals are identified as the origin of BPPV. There, calcium carbonate particles (CaCO 3 ) called otoliths, which are normally affixed to the canal walls, are detached by the aforementioned head or body changes. This extra mass floating in the endolymph causes an abnormal movement of the cupula, since these particles displace more volume of endolymph than usual. The brain misinterprets this displacement and sends erroneous information to the eyes, provoking a characteristic ocular nystagmus and the subsequent vertigo.
Dix and Hallpike [3] were pioneers in developing maneuvers which led to the detection of BPPV. These maneuvers consist of a series of consecutive head turns that trigger ocular responses in the patient, on the basis of which clinicians can determine whether a patient suffers from BPPV or not. Rabbit [4] developed a model which calculates the volume of endolymph displaced when the Dix and Hallpike maneuver is put into practice. In this model, a mathematical approximation consisting of a curve crossing the different angular positions was used. This is a theoretical model never validated with real data as far as the authors know. In this paper, the model is particularized to a specific real situation and an experimental plan is produced. In our case, a maneuver composed of two consecutive turns has been developed: a positive pitch (turning the head back) from the upright position, followed by a positive roll (turning the head right), as these are the most common two head movements that trigger the above-mentioned nystagmus. In order to Figure 1 Semicircular canals and ampulla. Layout of the semicircular canals and ampula, both located in the inner ear. http://www.tbiomed.com/content/10/1/21 reflect a more realistic situation, a cubic interpolation between points has been carried out. Another advantage of carrying out the cubic interpolation with respect to Rabbit's model, is to obtain an analytical expression of the curve. This will permit us to design statistical experiments aimed at deriving an optimal estimation of the unknown parameters of the model.

Optimal experimental design
One of the main aims of Statistics consists of modeling the behaviour of any stochastic process or particular system. In order to carry this out, a mathematical model is always needed to explain the results obtained once the process has occurred, and predict accurately the future behaviour of that process or system. The mathematical models used in real situations depend on unknown parameters. In order to describe the way in which the results are expected to vary, we need to estimate these parameters with optimal accuracy. For that purpose, the Design of Experiments is used to help us design how to change the inputs of processes in order to observe and identify the reasons of the changes observed in the response y, which is usually expressed as follows, where z represents the set of observations to collect and θ the unknown parameters of the model. The error ε follows a Gaussian distribution with mean 0 and constant variance σ 2 , that is E(y) = η(z, θ) and var(y) = σ 2 .
The problem of determining which set of observations to collect is what will define the design. It is common to say that the input variables are controlled by the researcher, while the unknown parameters are determined by nature. Optimal Design of Experiments theory allows us to find the best design in the sense of obtaining an optimal estimate of the parameters of the model. Next, basic concepts of this theory will be briefly presented as well as the two main criteria to obtain this optimal estimation. Following that, we will explain how the model for the maneuver has been constructed and finally, optimal designs for this model are calculated both for discrete and continuous design space.
An exact experimental design of size N consists of a set of N observations collected at points z 1 , . . . , z N , in a given compact space χ. Some of these N points may be repeated, meaning that several observations are taken at the same value of z. This number of observations is usually predetermined by experimental cost constraints. A convenient way to understand designs is to treat them as a collection of different points of χ, together with the proportion of the N observations to be allocated at each different point. This suggests the idea of extending the definition of a design to any probability measure ξ on χ (approximate design), that is The collection of weights, p i = ξ(z i ), provide a probability measure on χ supported on the points z 1 , . . . , z k . Thus, an experiment will be replicated about Np i times on value z i . Kiefer [5] pioneered this approach, and its many advantages are well documented in design monographs, see Silvey [6] for example. This approach has been applied to optimal http://www.tbiomed.com/content/10/1/21 treatment allocation [7], optimal estimation of kinetic parameters of the Michaelis-Menten model [8] and the Arrhenius equation [9]. In what follows, the approximate design approach is adopted without loss of generality, restricting the attention to designs with a finite set of support points. For convenience, the design is described using a tworow matrix, with the support points displayed in the first row and their corresponding proportions of observations in the second row (2).
Let θ T be the unknown parameter vector, and let be evaluated at the nominal value of the parameter θ (0) . This nominal value usually represents the best guess for the parameters θ at the beginning of the experiment and it is necessary to proceed as one would with a non-linear model [10]. Assuming model errors follow a normal distribution and apart from an unimportant multiplicative constant, the Information Matrix of a design ξ is given by where p i = ξ(z i ) is the proportion of observations to be taken at point z i (see e.g. [10]). The covariance matrix of the least squares and maximum likelihood estimatorθ is asymptotically proportional to the inverse of this matrix [6]. The use of this matrix is very important when it comes to designing the experiment in an optimal way. The objective to be achieved is to find a design which gives the best estimation of the parameters (or linear functions of them), usually by using the least squares method or maximum-likelihood estimation method. Through what it is defined as criteria , we will be able to measure the accuracy of the design and to compare different designs of the same model.
The design criteria used in this work for estimating the model parameters are Doptimality and c-optimality [11]. The D-optimality criterion minimizes the volume of the confidence ellipsoid of the parameters and is given by where m is the number of parameters in the model. The D-optimal design will be that which minimizes the function D [ M(ξ , θ)]. The c-optimality criteria is used to estimate a linear combination of the parameters, say c T θ; it is the variance of this estimate which It is known that these criteria are all convex and nonincreasing functions of the designs and so, designs with small criterion values are desirable [6]. A design that minimizes one of these functions over all the designs on χ is called a -optimal design, or more specifically, a D-or a c-optimal design, respectively. An advantage of working with approximate designs is that the optimality of a particular design can be checked mathematically. Since the criteria are convex, standard convex analysis arguments using directional derivatives when is differentiable [6], show that a design ξ * is -optimal, if and only if, it satisfies (Equivalence Theorem) where (ξ * ) is [ M(ξ * , θ)] for short and ∇ (ξ * ) denotes the gradient of (ξ * ). The equality is reached at the support points of ξ * . When = D , that is, the D-optimality where m is the number of parameters. Therefore, the inequality (4) becomes is known as the generalized variance. This important theorem provides methods for constructing optimal designs [11,12].

Derivation of a model for the maneuver
The model to be derived will be used to predict the endolymph volume displacement in response to a maneuver composed by two consecutive turns: a positive pitch (turning the head back) from the upright position followed by a positive roll (turning the head right). The standard definitions are used for head rotations: pitch rotates about an axis out the ear (y-axis), roll rotates about an axis out the eye (x-axis) and yaw rotates about an axis pointing out the top of the head (z-axis). We are only concerned with the horizontal canal since usually the particles are located inside that area [13]. Therefore, considering Q( t, α, θ) to be the volume of endolymph displaced after the first turn, where t and α are the variables which represent duration and angle of each turn, respectively, our model is formulated as follows, The error ε follows a Gaussian distribution with mean 0 and constant variance σ 2 , that is The forces acting on the endolymph due to the viscous interaction with a free-floating particle inside the canal are the interaction drag force F between the particle and endolymph and the gravity force g acting on the particle. After considering the forces mentioned, as well as the simplifications corresponding to the fluid, the Navier-Stokes equations must be applied with the corresponding boundary and initial conditions [14]. Thus, the equations which relate the volume flow rate of endolymph to the pressure and inertial forces at time t are The unknown parameters, θ T = (θ 1 , θ 2 ), are related to the damping and stiffness of the canal, respectively. The right side of the first equation is represented by two forces. The inertial forcing due to the angular acceleration of the head-fixed system relative to the inertial frame (ground-fixed system) is where¨ (t) ∧ R(s) stands for the tangential acceleration,¨ (t) being the angular acceleration of the head relative to the ground-fixed inertial frame resolved into the head-fixed frame and R(s) the vector running from the head-fixed coordinate system's origin to the centerline of the canal. The parameterization of R(s) is made with respect to the arc length s (also called natural parameter). The head-fixed coordinate system was defined when the subject was in the upright position prior to movement of the head. The constant ρ stands for the density of the canal and l n is the length covered by the otolith. http://www.tbiomed.com/content/10/1/21 The second term, F n , is the result of the interaction drag forces due to the particle moving relative to the fluid, In this equation, g is the gravitational acceleration, n is the unit normal tangent vector to the canal centerline andξ is the velocity of the particle. The constants A s , N, A e , a, ρ s , ρ e and μ e , stand for frontal area of the particle, number of particles which are floating inside the canal, cross-sectional area of the canal, radius of the particle, density of the particle, density of the endolymph and endolymph viscosity, respectively.
The manner in which the maneuver has been carried out was established by the clinician. The angle of turn in each time t is represented by (t). The movements of the head determine the magnitude and direction of the vectors (t) and g. It is important to note that the model equations refer to the non-inertial system, therefore, the linear and angular acceleration must be resolved into this system. This is done by using where¨ I (t) is the angular acceleration referred to the inertial system (for example the clinician who makes the maneuvers) and M(t) a rotation matrix. Since these maneuvers consist of a pitch (y-axis) followed by a roll movement (x-axis), the vector I for the pitch is (0, (t), 0) T and for the roll it is ( (t), 0, 0) T .
The rotation matrix for the first turn is and for the second turn, Once boundary conditions for the angles are imposed, (t) is written as follows, In order to obtain the expressions for the forcing terms F i and F n , we assume that the geometry of the posterior canal is described by a circle of radius r. Therefore, and F n = a 2 N b 2 4 3 a(ρ s − ρ e ) −r sin α sin 2 (t)¨ 2 (t) − g cos(ξ/r) + 6 μξ a + a N 6 μ e π b 2 (b 2 − a 2 )Q , http://www.tbiomed.com/content/10/1/21 where A s = πa 2 , A e = πb 2 , b stands for the radius of the cross-sectional area of the canal andξ = 0.02 cm/s [2]. After imposing the initial condition, the solution of equations (6) is written as

Optimal experimental designs
Optimal designs for the model given in (5) with nominal values of the parameters θ (0) 1 = 0.85 and θ (0) 2 = 0.2 (obtained from [15]) are calculated. It is assumed that six particles are detached from the canal wall [16]. On the other hand, it is known that A s /A e ≈ 10 −4 and that, approximately, the values of r and l n are 0.1 cm and 0.5 cm, respectively [17]. The values used for the constants are listed in Table 1. Figure 2 shows the plot of the function given by the above mentioned Q( t, α, θ). This is plotted as a function of t and α for the given nominal values.
The Information Matrix is expressed as where The two components are written as and

D-optimal design for a discrete design space
D-optimal designs for model (5)  Taking this finite set into consideration and maximizing the determinant of the Information Matrix, the optimal design is ξ * 1 = ( 1.5 , π/6 ) ( 0.5 , π/4 ) ( 0.5 , π/2 ) 0.47 0.16 0.37 . Table 2 shows that all the values of the generalized variance are smaller than the number of parameters, and therefore the design obtained satisfies numerically the Equivalence Theorem. The design indicates the duration and angle of the maneuvers to be carried out, but not how many observations the sample will have. That will be determined by the researcher through other means (for instance, the budget). The weights will give us the proportion of different maneuvers to be performed. For instance, if the sample contains N = 100 observations, the clinician will have to repeat the maneuver consisting of a positive 90 • pitch from the upright position followed by a positive 90 • roll, both with a duration of half a second, 37 times.

Sensitivity analysis of the D-optimal design
The goodness of any design ξ is measured by its efficiency. We will see how efficient the D-optimal design ξ * 1 is by checking the value of where ξ * θ is the D-optimal design calculated for a possible real value of the parameter θ. This efficiency shows how robust the design is with respect to the true (unknown) value of the parameters. To check the robustness of a design, we want to check how the quality of the estimation would be affected by a wrong choice of the nominal value. The efficiency can sometimes be multiplied by 100 and be reported in percentage terms. If, for instance, we take as nominal value θ (0) = 0.5, the true value being θ = 0.6 and the efficiency being 50%, then the design ξ * θ (0) would need to double the total number of observations to perform as well as the optimal design calculated with the true value θ = 0.6. Thus, our design would not be very robust. Table 3 illustrates the sensitivity of the D-optimal design with respect to the choice of the parameters θ 1 (horizontal values) and θ 2 (vertical values), i.e, how robust the design is with respect to the true values of those parameters. As can be observed, small variations of the nominal values does not affect the quality of the estimation much.

D-optimal design for a continuous design space
If instead of the finite set χ = {π/6 , π/4 , π/3 , 5π/12 , π/2} × {0.5 , 1 , 1.5}, we would have chosen a continuous design space χ =[ 0.5, 1.5] ×[ π/6, π/2], the Doptimal design obtained would have been As can be observed, the results would have been very similar to ξ 1 , having an efficiency of around 95%. But in practice, the clinician cannot carry out any turns that are smaller than 15°or shorter than half a second. Figure 3 shows how the generalized variance function for a continuous design space satisfies the Equivalence Theorem, that is, for all the values of t and α within the space design, the values of the generalized variance must be lower or equal to the number of parameters to estimate. At the points of the design, the value of the function must be equal to the number of parameters.

c-optimal design for a discrete design space
If we are interested in estimating a linear combination of the parameters, say c T θ, then we use the c-optimality criterion. An elegant way for finding an optimal design that estimates The generalized variance for the design ξ * 2 verifies the Equivalence Theorem. http://www.tbiomed.com/content/10/1/21 Figure 3 Plot of the generalized variance. Plot of the generalized variance for the design obtained. This design will be D-optimal to estimate both parameters θ 1 and θ 2 , if and only if, it verifies the Equivalence Theorem. That is, for all the values of t and α within the space design, the values of the generalized variance must be lower than or equal to than the number of parameters to estimate. In the points of the design, the value of the function must be equal to the number of parameters.
a linear combination of the parameters was given by Elfving [18]. This method is nicely illustrated and explained, e.g. by Chernoff [19], Kitsos [20] or Wiley [10]. For a given regression problem with regression function f ( t, α, θ), the method first defines the Elfving's set given by set G, the convex hull of {f ( t, α, θ (0) ) ∪ −f ( t, α, θ (0) )}, θ (0) being a nominal value for the parameter θ. This means that the set G is the smallest convex set containing {f ( t, α, θ (0) ) ∪ −f ( t, α, θ (0) )}. The point of intersection of the straight line defined by the vector c with the boundary of the Elfving's set determines the c-optimal design, ξ * , as a convex combination of the vertices of G. These vertices provide the support points of the optimal design. The weights in the convex combinations are the weights of the optimal design. Furthermore, where c * is the vector defined by the cut point of the straight line defined by c with the boundary of G.
On the other hand, to estimate the parameter θ 2 , let c T = (0, 1). The straight line defined by this vector cuts the boundary of the Elfving's set at point (0, y * ), and it defines the c θ 2 -optimal design to estimate the parameter θ 2 . This point is expressed as a convex combination of points f ( t 1 , α 1 ) and −f ( t 4 , α 4 ) = −f (1.5, π/6) with weights p 1 = 0.3 and p 2 = 0.7, respectively. Therefore, the optimal design is One important aspect is to check the efficiency of the D-optimal design (ξ * D ) with respect to the c-optimal design (ξ * c ). The formula of the efficiency, provides a way to see how good the D-optimal design is at estimating each of the parameters. In the case of θ 1 , the D-optimal design ξ * 1 and c-optimal design ξ * 3 are compared and the efficiency is around 75%, while for θ 2 , ξ * 1 and ξ * 4 are compared, having an efficiency of around 85%. These results show the D-optimal design is more efficient for estimating θ 2 than for estimating θ 1 , that is, with this design, the test power for testing {H 0 : θ 2 = 0} will be greater that the test power for {H 0 : θ 1 = 0}.

c-optimal design for a continuous design space
For a continuous design space, the Elfving's set obtained is shown in Figure 5 by plotting the convex hull of the surface {f ( t, α, θ (0) ) ∪ −f ( t, α, θ (0) )}. In this case t ∈ [ 0.5, 1.5] ×[ π/6, π/2]. With the help of the method proposed by López-Fidalgo and Rodríguez-Díaz [21], we will calculate the c-optimal design to estimate the parameters θ 1 , θ 2 and an example of a linear combination of them, θ 1 + θ 2 . The efficiency of the D-optimal design for estimating the parameter θ 1 is around 75%, while for θ 2 , it is around 90%. Finally, to obtain the c-optimal design for the estimation of θ 1 + θ 2 , vector c T = (1, 1) is considered. In this case, the straight line defines the point c * = (0.084, 0.084) = −f (1.3, 0.6), so the c-optimal design is supported at one point: As we can observe, the results obtained are quite similar to the case concerning the discrete design space, except for the estimation of θ 1 + θ 2 , where for the continuous case the c-optimal design is only supported at one point, although in both cases the angle of turn is similar. If the number of support points in a c-optimal design is less than the number of parameters, this design allows the computation of the maximum likelihood estimate of this linear combination. But, in this case, not all the parameters are identifiable individually, that is, some of them cannot be estimated.

Discussion
The present biomechanical model was derived to provide a quantitative basis for the detection of BPPV. This model is based on a maneuver consisting of two consecutive head turns. These are the most common head movements leading to vertigo symptoms. We would like to remark that although the model can only be applied for this specific maneuver, it could be extended to other types. The experiment, that is, the duration and angle of the head movements to be applied to the patients should be based on the design provided http://www.tbiomed.com/content/10/1/21 by the D-optimal design since it helps us estimate the parameters simultaneously, minimizing the confidence ellipsoid. The c-optimal design is used either for estimating linear combinations of the parameters, or for estimating the parameters separately. But in this case, it also provides valuable assistance to check how efficient the D-optimal design is for the estimation of each of the parameters. This is an interesting check of the sensitivity since a D-optimal design could be quite efficient for estimating a particular parameter but quite inefficient for estimating another one.
The covariance matrix of the estimates is asymptotically proportional to the inverse of the Fisher Information Matrix. Theoretically, the application of the maneuvers which are specified in the design, along with their corresponding proportions, will assure that an objective function of the covariance matrix of the estimators θ = varθ 1 cov (θ 1 ,θ 2 ) cov (θ 1 ,θ 2 ) varθ 2 ∝ M −1 (ξ , θ), will be minimized. Symbol ∝ stands for "asymptotically proportional" in this case. Finally, we would like to point out that, as far as the authors know, the models found in the published works describing this sort of maneuvers have not been validated with data yet. The clinicians hold that the extra volume of endolymph displaced by the otoliths are directly related to the eye movements provoked in the patient under vertigo symptoms. Therefore, in some way, to validate the model, response y should be measured through some variable related to eye movement.