extfieldsofstudy:
Computer Science
Physics
Mathematics

text:
Curvature-Exploiting Acceleration of Elastic Net Computations

This paper introduces an efficient second-order method for solving the elastic net problem. Its key innovation is a computationally efficient technique for injecting curvature information in the optimization process which admits a strong theoretical performance guarantee. In particular, we show improved run time over popular first-order methods and quantify the speed-up in terms of statistical measures of the data matrix. The improved time complexity is the result of an extensive exploitation of the problem structure and a careful combination of second-order information, variance reduction techniques, and momentum acceleration. Beside theoretical speed-up, experimental results demonstrate great practical performance benefits of curvature information, especially for ill-conditioned data sets.

Introduction
Lasso, ridge and elastic net regression are fundamental problems in statistics and machine learning, with countless applications in science and engineering [33]. Elastic net regression amounts to solving the following convex optimization problem minimize x∈Ê d

2n
Ax for given data matrices A ∈ Ê n×d and b ∈ Ê n and regularization parameters γ 1 and γ 2 .
Setting γ 1 = 0 results in ridge regression, γ 2 = 0 yields lasso and letting γ 1 = γ 2 = 0 reduces the problem to the classical least-squares. Lasso promotes sparsity of the optimal solution, which sometimes helps to improve interpretability of the results. Adding the additional l 2 -regularizer helps to improve the performance when features are highly correlated [28,33].
The convergence rates of iterative methods for solving (1) are typically governed by the condition number of the Hessian matrix of the ridge loss, C + γ 2 I, where C = 1 n A ⊤ A is the sample correlation matrix. Real-world data sets often have few dominant features, while the other features are highly correlated with the stronger ones [13,28]. This translates to a rapidly decaying spectrum of C. In this paper, we demonstrate how this property can be exploited to reduce the effect of ill-conditioning and to design faster algorithms for solving the elastic net regression problem (1).

Related work
Over the past few years, there has been a great attention to developing efficient optimization algorithms for minimizing composite objective functions where f (x) = 1 n n i=1 f i (x) is a finite sum of smooth and convex component functions f i (x), and h (x) is a possibly non-smooth convex regularizer. In machine learning applications, the function f typically models the empirical data loss and the regularizer h is used to promote desired properties of a solution. For example, the elastic net objective can fit to this form with f i (x) = 1 2 a ⊤ i x − b i 2 + γ 2 x 2 2 , and h(x) = γ 1 x 1 .

First-Order methods
Standard deterministic first-order methods for solving (2), such as proximal gradient descent, enjoy linear convergence for strongly convex objective functions and are able to find an ǫapproximate solution in time O dnκ log 1 ǫ , where κ is the condition number of f . This runtime can be improved to O dn √ κ log 1 ǫ if it is combined with Nesterov acceleration [4,23]. However, the main drawback of these methods is that they need to access the whole data set in every iteration, which is too costly in many machine learning tasks.
For large-scale problems, methods based on stochastic gradients have become the standard choice for solving (2). Many linearly convergent proximal methods such as, SAGA [9] and Prox-SVRG [31], have been introduced and shown to outperform standard first-order methods under certain regularity assumptions. These methods improve the time complexity to O d (n +κ) log 1 ǫ , whereκ is a condition number satisfyingκ ≥ κ. When the component functions do not vary substantially in smoothness,κ ≈ κ, and this complexity is far better than those of deterministic methods above. By exploiting Nesterov momentum in different ways (see, e.g., [2,8,11,19]), one can improve the complexity to O d n + √ nκ log 1 ǫ , which is also optimal for this class of problems [30].

Second-order methods
Second-order methods are known to have superior performance compared to their first-order counterparts both in theory and practice, especially when the problem at hand is highly nonlinear and/or ill-conditioned. However, such methods often have very high computational cost per iteration. Recently, there has been an intense effort to develop algorithms which use second-order information with a more reasonable computational burden (see, e.g., [1,6,10,21,26,31,32] and references therein). Those methods use techniques such as random sketching, matrix sampling, and iterative estimation to construct an approximate Hessian matrix. Local and global convergence guarantees have been derived under various assumptions. Although many experimental results have shown excellent performance of those methods on many machine learning tasks, current second-order methods for finite-sum optimization tend to have much higher time-complexities than their first-order counterparts (see [32] for a detailed comparison).
Apart from having high time complexities, none of the methods cited above have any guarantees in the composite setting since their analyses hinge on differentiability of the objective function. Instead, one has to rely on methods that build on proximal Newton updates (see, e.g., [12,18,20,25]). However, these methods still inherit the high update and storage costs of conventional second-order methods or require elaborate tuning of several parameters and stopping criteria depending on a phase transition which occurs in the algorithm.

Ridge regression
For the smooth ridge regression problem, the authors in [13] have developed a preconditioning method based on linear sketching which, when coupled with SVRG, yields a significant speed-up over stochastic first-order methods. This is a rare second-order method that has a comparable or even better time complexity than stochastic first-order methods. More precisely, it has a guaranteed running time of O(d(n + κ H ) log 1 ǫ ), where κ H is a new condition number that can be dramatically smaller thanκ, especially when the spectrum of C decays rapidly. When d ≪ n, the authors in [29] combine sub-sampled Newton methods with the mini-batch SVRG to obtain some further improvements.

Contributions
Recently, the work [3] shows that under some mild algorithmic assumptions, and if the dimension is sufficiently large, the iteration complexity of second-order methods for smooth finite-sum problems composed of quadratics is no better than first-order methods. Therefore, it is natural to ask whether one can develop a second-order method for solving the elastic net problem which has improved practical performance but still enjoys a strong worst-case time complexity like the stochastic first-order methods do? It should be emphasized that due to the non-smooth objective, achieving this goal is much more challenging than for ridge regression. The preconditioning approach in [13] is not applicable, and the current theoretical results for second-order methods are not likely to offer the desired running time.
In this paper, we provide a positive answer to this question. Our main contribution is the design and analysis of a simple second-order method for the elastic net problem which has a strong theoretical time complexity and superior practical performance. The convergence bound adapts to the problem structure and is governed by the spectrum and a statistical

Algorithm
Time complexity 2nd-order yes measure of the data matrix. These quantities often yield significantly stronger time complexity guarantees for practical datasets than those of stochastic first-order methods (see Table 1). To achieve this, we first leverage recent advances in randomized low-rank approximation to generate a simple, one-shot approximation of the Hessian matrix. We then exploit the composite and finite-sum structure of the problem to develop a variant of the ProxSVRG method that builds upon Nesterov's momentum acceleration and inexact computations of scaled proximal operators, which may be of independent interest. We provide a simple convergence proof based on an explicit Lyapunov function, thus avoiding the use of sophisticated stochastic estimate sequences.

Preliminaries and Notation
Vectors are indicated by bold lower-case letters, and matrices are denoted by bold uppercase letters. We denote the dot product between x and y by x, y = x ⊤ y, and the Euclidean norm of x by x 2 = x, x . For a symmetric positive definite matrix H, x, y H = x ⊤ Hy is the H-inner product of two vectors x and y and x H = x, x H is the Mahalanobis norm of x. We denote by λ i (A) the ith largest eigenvalue of A. Finally, λ i denotes the ith largest eigenvalue of the correlation matrix C.
In the paper, we shall frequently use the notions of strong convexity and smoothness in the H-norm, introduced in the next two assumptions.
Assumption 1 implies that ∇f is L-Lipschitz: As a consequence, we have the following bound: Assumption 2. The function f (x) is µ-strongly convex w.r.t the H-norm, i.e, there exists a positive constant µ such that holds for all x, y ∈ Ê n and λ ∈ [0, 1].
Assumption 2 is equivalent to the requirement that holds for all x, y ∈ Ê n . We will use both of these definitions of strong convexity in our proofs.
At the core of our method is the concept of scaled proximal mappings, defined as follows: Definition 1 (Scaled Proximal Mapping). For a convex function h and a symmetric positive definite matrix H, the scaled proximal mapping of h at x is The scaled proximal mappings generalize the conventional ones: However, while many conventional prox-mappings admit analytical solutions, this is almost never the case for scaled proximal mappings. This makes it hard to extend efficient firstorder proximal methods to second-order ones. Fortunately, scaled proximal mappings do share some key properties with the conventional ones. We collect a few of them in the following result: ). The following properties hold: 1. prox H h (x) exists and is unique for x ∈ dom h. 2. Let ∂h (x) be the subdifferential of h at x, then Finally, in our algorithm, it will be enough to solve (3) approximately in the following sense: Definition 2 (Inexact subproblem solutions). We say that x + ∈ Ê d is an ǫ-optimal solution 3 Building Block 1: Randomized Low-Rank Approximation The computational cost of many Newton-type methods is dominated by the time required to compute the update direction d = H −1 g for some vector g ∈ Ê d and approximate Hessian H. A naive implementation using SVD would take O (nd 2 ) flops, which is prohibitive for large-scale data sets. A natural way to reduce this cost is to use truncated SVD. However, standard deterministic methods such as the power method and the Lanzcos method have run times that scale inversely with the gap between the eigenvalues of the input matrix. This gap can be arbitrarily small for practical data sets, thereby preventing us from obtaining the desired time complexity. In contrast, randomized sketching schemes usually admit gapfree run times [15]. However, unlike other methods, the block Lanczos method, detailed in Algorithm 1, admits both fast run times and strong guarantees on the errors between the true and the computed approximate singular vectors. This property turns out to be critical for deriving bounds on the condition number of the elastic net.
i be the SVD of A. Then, the following bounds hold with probability at least 9/10: Note that we only run Algorithm 1 once and ǫ ′ = 1/2 is sufficient in our work. Thus, the computational cost of this step is negligible, in theory and in practice.

Aproximating the Hessian
In this work, we consider the following approximate Hessian matrix of the ridge loss: Here, the first term is a natural rank r approximation of the true Hessian, while the second term is used to capture information in the subspace orthogonal to the column space of V r . The inverse of H in (6) admits the explicit expression so the evaluation of H −1 x has time complexity O (rd).

Bounding the Condition Number
We now turn our attention to studying how the approximate Hessian affects the relevant condition number of the elastic net problem. We first introduce a condition number that usually determines the iteration complexity of stochastic first-order methods under nonuniform sampling.
Definition 3. The average condition number of (1.1) is For the elastic net problem (1), the smooth part of the objective is the ridge loss .
Since we define smoothness and strong convexity of f i (x) in the H-norm, the relevant constants are For comparison, we also define the conventional condition numberκ, which characterizes the smoothness and strong convexity of f i (x) in the Euclidean norm. In this caseκ = and µ = λ d (C + γ 2 I) .
It will become apparent that κ H can be expressed in terms of a statistical measure of the ridge loss and that it may be significantly smaller thanκ. We start by introducing a statistical measure that has been widely used in the analysis of ridge regression (see, e.g., [16] and the references therein).

Definition 4 (Effective Dimension).
For a positive constant λ, the effective dimension of C is defined as The effective dimension generalizes the ordinary dimension and satisfies d λ ≤ d with equality if and only if λ = 0. It is typical that when C has a rapidly decaying spectrum, most of the λ i 's are dominated by λ, and hence d λ can be much smaller than d.
The following lemma bounds the eigenvalues of the matrix H −1/2 (C + γ 2 I) H −1/2 , which can be seen as the effective Hessian matrix.

Lemma 1 ([13]
). Invoking Algorithm 1 with data matrix 1 √ n A, target rank r, and target precision ǫ ′ = 1/2, it holds with probability at least 9/10 that Equipped with Lemma 1, we can now connect κ H with d λ using the following result.
Theorem 1. With probability at least 9/10, the following bound holds up to a multiplicative constant: Proof. See Appendix A.
compared toκ. If the spectrum of C decays rapidly, then the terms i>r λ i are negligible and the ratio is approximately i≤r λ i /(rλ r ). If the first eigenvalues are much larger than λ r , this ratio will be large. For example, for the australian data set [7], this ratio can be as large as 1.34 × 10 4 and 1.6 × 10 5 for r = 3 and r = 4, respectively. This indicates that it is possible to improve the iteration complexity of stochastic first-order methods if one can capitalize on the notions of strong convexity and smoothness w.r.t the H-norm in the optimization algorithm. Of course, this is only meaningful if there is an efficient way to inject curvature information into the optimization process without significantly increasing the computational cost. In the smooth case, i.e., γ 1 = 0, this task can be done by a preconditioning step [13]. However, this approach is not applicable for the elastic net, and we need to make use of another building block.

Building Block 2: Inexact Accelerated Scaled Proximal SVRG
In this section, we introduce an inexact scaled accelerated ProxSVRG algorithm for solving the generic finite-sum minimization problem in (2). We then characterize the convergence rate of the proposed algorithm.

Description of the Algorithm
To motivate our algorithm, we first recall the ProxSVRG method from [31]: For the sth outer iteration with the corresponding outer iteratex s , let x 0 =x s and for k = 0, 2, . . . , where i k is drawn randomly from {1, . . . , n} with probability p i k = L i k /(nL avg ). Since we are provided with an approximate Hessian matrix H, it is natural to use the following update: which can be seen as a proximal Newton step with the full gradient vector replaced by v k . Note that when h(·) is the ℓ 1 -penalty, ProxSVRG can evaluate (8) in time O(d), while evaluating (9) amounts to solving an optimization problem. It is thus is critical to keep the number of such evaluations small, which then translates into making a sufficient progress at each iteration. A natural way to achieve this goal is to reduce the variance of the noisy 10: end for 11:x s+1 ← x T 12: end for Output:x S gradient v k . This suggests to use large mini-batches, i.e., instead of using a single component function f i k , we use multiple ones to form: where B k ⊂ {1, . . . , n} is a set of indices with cardinality |B k | = b. It is easy to verify that v k is an unbiased estimate of ∇f (x k ). Notice that naively increasing the batch size makes the algorithm increasingly similar to its deterministic counterpart, hence inheriting a high-time complexity. This makes it hard to retain the runtime of ProxSVRG in the presence of 2nd-order information.
In the absence of second-order information and under the assumption that the proximal step is computed exactly, the work [24] introduced a method called AccProxSVRG that enjoys the same time complexity as ProxSVRG but allows for much larger mini-batch sizes. In fact, it can tolerate a mini-batch of size O √κ thanks to the use of Nesterov momentum. This indicates that an appropriate use of Nesterov momentum in our algorithm could allow for the larger mini-batches required to balance the computational cost of using scaled proximal mappings. The improved iteration complexity of the scaled proximal mappings will then give an overall acceleration in terms of wall-clock time. As discussed in [2], the momentum mechanism in AccProxSVRG fails to accelerate ProxSVRG unlessκ ≥ n 2 . In contrast, as we will see, our algorithm will be able to accelerate the convergence also in these scenarios. In summary, our algorithm is designed to run in an inner-outer fashion as ProxSVRG with large mini-batch sizes and Nesterov momentum to compensate for the increased computational cost of subproblems. The overall procedure is summarized in Algorithm 2.

Convergence Argument
In this subsection, we will show that as long as the errors in evaluating the scaled proximal mappings are controlled in an appropriate way, the iterates generated by the outer loop of Algorithm 2 converge linearly in expectation to the optimal solution. Recall that in Step 7 of Algorithm 2, we want to find an ǫ k -optimal solution in the sense of (5) to the following problem: The next lemma quantifies the progress made by one inner iteration of the algorithm. Our proof builds on a Lyapunov argument using a Lyapunov function on the form: Proof. See Appendix C.
Remark 1. Our proof is direct and based on natural Lyapunov functions, thereby avoiding the use of stochastic estimate subsequences as in [24] which is already very complicated even when the subproblems are solved exactly and H = I. We stress that the result in Lemma 2 also holds for smaller mini-batch sizes, namely b ∈ 1, . . . , O( L avg /µ) , provided that the step size η is reduced accordingly. In favor of a simple proof, we only report the large mini-batch result here.
Equipped with Lemma 2, we can now characterize the progress made by one outer iteration of Algorithm 2.
Proof. See Appendix D.
Remark 2. The theorem indicates that if the errors in solving the subproblems are controlled appropriately, the outer iterates generated by Algorithm 2 converge linearly in expectation to the optimal solution. Since V 0 depends on x ⋆ , it is difficult to provide a general closed-form expression for the target precisions ǫ k . However, we will show below that with a certain policy for selecting the initial point, it is sufficient to run the solver a constant number of iterations independently of x ⋆ . We stress that the results in this section are valid for minimizing general convex composite functions (2) and not limited to the elastic net problem.

Warm-Start
The overall time complexity of Algorithm 2 depends strongly on our ability to solve problem (10) in a reasonable computational time. If one naively starts the solver at a random point, it may take many iterations to meet the target precision. Thus, it is necessary to have a well-designed warm-start procedure for initializing the solver. Intuitively, the current iterate x k can be a reasonable starting point since the next iterate x k+1 should not be too far away from x k . However, in order to achieve a strong theoretical running time, we use a rather different scheme inspired by [19]. Let us first define the vector u k = y k − ηH −1 v k for k ∈ {0, 1, . . . , T − 1} and the function Then, the kth subproblem seeks for x k+1 such that where x ⋆ k+1 is the exact solution. We consider the initialization policy which can be seen as one step of the proximal gradient method applied to p (z, u k−1 ) starting at the current x k with step size γ.
The following proposition characterizes the difference in objective realized by z 0 and x ⋆ k+1 .
The proposition, together with (13), implies that it suffices to find x k+1 such that This is significant since one only needs to reduce the residual error by a constant factor independent of the target precision. Note also that the condition number κ sub is much smaller than κ ≈ λ 1 (H)/(λ d + γ 2 ), and computing the gradient of the smooth part of p (z, u) only takes time O (rd) instead of O (nd) as in the original problem. Those properties imply that the subproblems can be solved efficiently by iterative methods, where only a small (and known) constant number of iterations is needed. The next section develops the final details of convergence proof.

Global Time Complexity
We start with the time complexity of Algorithm 2. Let T (α) be the number of gradient evaluations that a subproblem solver takes to reduce the residual error by a factor α. Then, by Proposition 2, one can find an ǫ k -optimal solution to the kth subproblem by at most T (κ sub / (1 − ρ)) gradient evaluations, where one gradient evaluation is equivalent to O (d) flops. Consider the same setting of Theorem 2 and suppose that the subproblems are initialized by (14). Then, the time complexity of Algorithm 2 is given by: where the first summand is due to the full gradient evaluation at each outer loop; the second one comes from the fact that one needs O( √ κ H ) inner iterations, each of which uses a minibatch of size O( √ κ H ); and the third one is the result of O( √ κ H ) inner iterations, each of which solves a subproblem that needs T (κ sub / (1 − ρ)) gradient evaluations. We can now put things together and state our main result.
Proposition 3. Suppose that the approximate Hessian matrix H is given by (6) and that Algorithm 2 is invoked with f (x) = 1 2n Ax − b 2 2 + γ 2 2 x 2 2 and h (x) = x 1 . Assume further that the subproblems are solved by the accelerated proximal gradient descent method [4,23]. Our method can find an ǫ-optimal solution in time Proof. The task reduces to evaluating the term T (κ sub / (1 − ρ)) in (16). Recall that the iteration complexity of the accelerated proximal gradient descent method for minimizing the function F (x) = f (x) + h(x), where f is a smooth and strongly convex function and h is a possibly non-smooth convex regularizer, initiliazed at x 0 is given by where κ is the condition number. By invoking the about result with F (x) = p (x, u k ), k+1 , x 0 = z 0 , κ = κ sub , and ǫ is the right-hand side of (15), it follows that the number of iterations for each subproblem can be bounded by In addition, each iteration takes time O (rd) to compute the gradient implying the time complexity Finally, since κ H ≫ κ sub , the proof is complete.
We can easily recognize that this time complexity has the same form as the stochastic first-order methods discussed in Section 1.2.1 with the condition numberκ replaced by κ H . It has been shown in Theorem 1 that κ H can be much smaller thanκ, especially, when C has a rapidly decaying spectrum. We emphasize that the expression in (17) is available for free to us after having approximated the Hessian matrix. Hence no tuning is required to set ǫ k for solving the subproblems.

Experimental Results
In this section, we perform numerical experiments to verify the efficacy of the proposed method on real world data sets [7,14]. We compare our method with several well-known first-order methods: FISTA [4] with optimal step-size; Prox-SVRG [31] with epoch length 2n/b as suggested by the authors; Katyusha1 [2] with epoch length 2n/b, Katyusha momentum τ 2 = 0.5/b as suggested by the author; and our method with epoch length 2n/b. Since Katyusha1 can use a mini-batch of size √ n without slowing down the convergence, we set b = √ n for all methods except FISTA. Finally, to make a fair comparison, for each algorithm above, we tune only the step size, from the set η × {10 k , 2 × 10 k , 5 × 10 k |k ∈ {0, ±1, ±2}}, where η is the theoretical step size, and report the one having smallest objective value. Other hyper-parameters are set to their theory-predicted values. For Katyusha1, we also compute the largest eigenvalue of the Hessian matrix in order to set its parameter τ 1 . All methods are initialized at 0 and run for up to 100 epochs (passes through the full data set). For the subproblems in Algorithm 2, we just simply run FISTA with √ κ sub log κ sub iterations as discussed in the previous section, without any further tunning steps. The value of r is chosen as a small fraction of d so that the preprocessing time of Algorithm 1 is negligible. Note that the available spectrum of C after running Algorithm 1 also provides an insightful way to choose r. Figure 2 shows the suboptimality in objective versus the number of epochs for different algorithms solving the elastic net problem. We can see that our method systematically outperforms the others in all settings, and that there is a clear correspondence between the spectrum of C in Fig. 6    the first-order methods make almost no progress in the first 100 epochs, while our method can find a high-accuracy solution within tens of epochs, demonstrating a great benefit of second-order information. On the other hand, for a well-conditioned data set that does not exhibit high curvature such as real-sim, ProxSVRG is comparable to our method and even outperforms Katyusha1. This agrees with the theoretical time complexities summarized in Table 1. Finally, we demonstrate a hard instance for low-rank approximation methods via the cina0 data set, where C has a very large condition number and slowly decaying dominant eigenvalues. In this case, low-rank approximation methods with a small approximate rank r may not be able to capture sufficient curvature information. We can see that even in this case curvature information helps to avoid the stagnation experienced by FISTA and ProxSVRG.

Conclusions
We have proposed and analyzed a novel second-order method for solving the elastic net problem. By carefully exploiting the problem structure, we demonstrated that it is possible to deal with the non-smooth objective and to efficiently inject curvature information into the optimization process without (significantly) increasing the computational cost per iteration. The combination of second-order information, fast iterative solvers, and a well-designed warm-start procedure results in a significant improvement of the total runtime complexity over popular first-order methods. An interesting direction for future research would be to go beyond the quadratic loss. We believe that the techniques developed in this work can be extended to more general settings, especially when the smooth part of the objective function is self-concordant.
[17] Chonghai Hu, James T. Kwok, and Weike Pan. Accelerated gradient methods for stochastic optimization and online learning. In Advances in Neural Information Processing Systems, pages 781-789, Vancouver, Canada, Dec 2009.

A Proof of Theorem 1
First, by the definition of L i , it holds that where the last inequality follows since The term a ⊤ i H −1 a i can be bounded as It can be verified that where we have used Lemma 1 in the last step. It follows that which implies that On the other hand, it was shown in [13, Theorem 4] that Finally, recall that µ = λ d H −1/2 (C + γ 2 I) H −1/2 . By Lemma 1, we readily have µ ≥ c 1 γ 2 for some constant c 1 , which concludes the proof.

B Some Useful Auxiliary Results
To facilitate the analysis, we collect some useful inequalities regarding the Mahalanobis norm that are used in the subsequent proofs.

C Proof of Lemma 2
Before proving the lemma, we rewrite the sequences {x k , y k , z k } k≥0 generated by Algorithm 2 in the following recurrence form: We also recall the following well-known three-point identity: which holds for any symmetric matrix H ≻ 0 and any vectors a, b, c ∈ Ê d .
We are now ready to prove Lemma 2. From the definition of the Lyapunov function in (11), we have that where the equality (a) follows from (19), and (b) follows from(18c). Using the identity (19) again for the last term in (20), we obtain where the last equality follows from (18b). By noting that where we have used (18a) and the fact that With these observations, we are now ready to bound F (x k+1 ). In particular, by invoking Lemma 3 with y = y t , (25), it follows that Thus, combining (27)- (30) and using the fact that Invoking Lemma 4 with u k = (1 − τ ) x k + τ x ⋆ and then applying Lemma 5, we obtain Therefore, by taking the expectation on both sides of (31) and using (32), we obtain By choosing η = 1 Lavg , τ = µ 2Lavg , and b ≥ 60 Lavg µ , it is readily verified that the second and third terms on the right-hand side of (33) become nonpositive, which concludes the proof.

D Proof of Theorem 2
To begin with, let e t = E {F (x t ) − F (x ⋆ )} andẽ s = E {F (x s ) − F (x ⋆ )}, then by applying inequality (12) recursively, we obtain It can be verified that Thus, from (34) and the facts that τ = µ 2Lavg , ξ t H −1 ≤ 2ǫ t /L avg , it holds that where By the definition of V k , we have µ Multiplying both sides of the above inequality by (1 − τ ) −k yields H , then we can write the previous inequality as where we have separated the term α 0 u 0 from the sum and added the positive term α k u k in the last step. Note that α 0 = 2 can be further bounded by where c 1 = 1 + √ 2 1−τ . It is readily verified that {B k } k is an increasing sequence and that {S k 2 µ (c 1 V 0 + B k )} k is an increasing sequence satisfying S 0 > u 2 0 . Therefore, by invoking Lemma 7 with u t , α k , and S k , we have for any k ≥ 1 that where we have used √ a + b ≤ √ a + √ b for any a, b ≥ 0. Thus, for any t ∈ {0, 1, . . . , k}, it holds that Having upper bounds of E z t − x ⋆ H , we can now substitute them into (35) to get where the equality follows from the definition of α t , and in the last step, we have added a positive term α k in the first sum as well as added and subtracted α 0 in the second sum. Note that the function √ 1 + x − x is decreasing on the interval [0, ∞), it follows that We thus have where in the last step, we have applied the inequality (a + b) 2 ≤ (1 + β) a 2 + 1 + 1 β b 2 twice with β = 9 and β = 1, respectively. Using the definitions of B k and α t , the above inequality can be rewritten as With our choice of ǫ t , it follows that where the third equality follows since by definition, e k−1 ≤ E V k−1 ≤ δ k−1 , and sinceẽ s = e 0 ≤ V 0 , and the last inequality follows since ρ = 0.9τ . By applying the above inequality recursively, we obtain where c 3 = c 2 + 2/9. Therefore, when T ≥ 4 3ρ log (27c 3 ), it holds that It follows that where the first inequality follows from the triangle inequality and the last one follows from Property 1. Therefore, it holds that Denote the filtration F k by F k = {f 0 , f 1 , . . . , f k }, where f i denotes all the randomness incurring at time i for all i ∈ {1, . . . , k}. Note that the triple (x k , y k , z k ) depends on F k−1 , but not on f k , we thus have Taking the expectation on both sides of (43) completes the proof.
Lemma 5 (Bounding Variance). Assume that the indices in the mini-batch B k are sampled independently from {1, . . . , n} with probabilities P i = L i nLavg , then conditioned on x k and y k , it holds that Proof. We follow the original proofs in [24,31] with a few modifications to take the Mahalanobis norm into account. For any i ∈ {1, . . . , n}, consider the function φ i defined by Then, φ i (x ⋆ ) = min x φ i (x) since ∇φ i (x ⋆ ) = 0 and φ i is convex. It can be checked that ∇φ i is L i -Lipschitz, we thus have which is equivalent to Multiplying both sides of the previous inequality by Lavg L i and averaging from 1, . . . , n gives


Secure Client and Server Geolocation Over the Internet

In this article, we provide a summary of recent efforts towards achieving Internet geolocation securely, \ie without allowing the entity being geolocated to cheat about its own geographic location. Cheating motivations arise from many factors, including impersonation (in the case locations are used to reinforce authentication), and gaining location-dependent benefits. In particular, we provide a technical overview of Client Presence Verification (CPV) and Server Location Verification (SLV)---two recently proposed techniques designed to verify the geographic locations of clients and servers in realtime over the Internet. Each technique addresses a wide range of adversarial tactics to manipulate geolocation, including the use of IP-hiding technologies like VPNs and anonymizers, as we now explain.

I. INTRODUCTION
Internet Geolocation is the process of determining the geographic location of an Internet-connected device. Secure geolocating of a web client (a client visiting a website) is useful for location-aware authentication, location-aware access control, location-based online voting, location-based social networking, and fraud reduction.
From the client's perspective, geolocating the remote server can provide higher assurance to the server's identity, and justify conducting certain sensitive transactions, e.g., those requiring certain privacy measures or requiring data sovereignty [1]. Independent of server and client geolocation, geolocating network intermediate systems (e.g., routers) can also be useful for purposes such as monitoring [2] and network mapping [3].
This article provides a technical overview on two recent approaches for secure location verification of clients and servers over the Internet, respectively CPV and SLV. Both techniques are based on network measurements, where delays are measured from trusted network nodes dubbed verifiers, and are analyzed in realtime to verify physical presence inside a prescribed geographic region. We explain the threat model of both techniques, how they mitigate against known adversarial tactics, how they adapt to various network dynamics, and what distinguishes them from other geolocation approaches.

II. GEOLOCATION BACKGROUND
This section provides a background on how Internet geolocation works, and the security limitations associated with common practices. Many academic geolocation methods have A version of this paper appeared in the USENIX ;login: magazine, Spring 2018 issue. This is the authors' copy for personal use. been proposed, but there have been very limited deployment in practice. As of this writing, most of the geolocation conducted in practice relies on the clients' IP address or GPS coordinates of hand-held devices, as explained below.

A. Geolocation in Practice
There are several methods for device geolocation over the Internet. If the device belongs to a user that is acting as a web client (i.e., visiting a website), the Geolocation API is a W3C standard that enables browsers to obtain location information of the device they are running on, and communicate it to a webserver. Servers request location coordinates using Javascript as follows: if(navigator.geolocation) { navigator.geolocation.getCurrentPosition(success, error , geoOptions); } else { console.log("Geolocation is not supported on this browser."); } The geolocation methods a browser uses are left to the browser vendor's discretion. Most major browsers rely on the following in varying orders (i.e., when one fails, the next is tried): GPS, WiFi Positioning System (WPS), IP address-based location look-ups, or cell-tower triangulation of mobile devices. The location of an IP address can be obtained from, e.g., publicly available routing information or public registries, such as whois. Many IP location service providers (commercial and free) maintain look-up tables to instantly map IP addresses to locations. Such static tablulation methods may take long times to reflect changes or IP address re-assignments, which occur quite often for client geolocation to be up to date (studies were conducted to confirm this [4]). IP address based geolocation can however be reliable for benign server geolocation. Flagfox is an example Firefox extension that visually indicates a flag of the country corresponding to the IP address resolution of the URL (see Fig 1). From a security perspective, none of the above techniques is resilient to adversarial manipulation. When the geolocation API is in use, the server normally makes no effort in geolocating the client device; it rather trusts the browser-communicated coordinates, which can easily be forged on the fly before being sent to the server. Example Firefox extensions that enable forgery include Fake Location (see Fig 2) and Location Guard; both enable a user to specify where in the world they would like to appear to be at. If the server relies on tabulation methods to geolocate the client (instead of asking the browser for its coordinates), the common practice of clients hiding their own IP addresses behind proxies and anonymizers comes into play.

B. Geolocation in the Literature
A wide set of techniques can be used, mostly for a server to geolocate clients [5]. These enable a server to infer a client's geographic location from, e.g., hints obtained from browsergenerated HTTP headers (such as preferred language or time zone). Locations can also be obtained through crowd-sourcing, e.g., by interpolating a device's location from its proximity to nearby devices (phones) with known GPS locations.
Another class of Internet geolocation approaches is based on network measurements. Similar to GPS triangulations that are based on the delays between the receiver and satellites, measurement-based techniques also aim to locate devices (clients or server) by estimating their distance from landmarks in the network with known locations. These landmarks measure network delays from themselves to the device, typically identified by its IP address, and map these delays to geographic distances. The accuracy of such mapping however is not anywhere near that of mapping satellite delays to distances, and is thus the primary source of inaccuracies in such techniques. Still, measurement-based geolocation is generally considered more accurate than methods like tabulation-based geolocation.
From the security point of view, although most of the above methods are positioned as resilient to evasion, examination has shown otherwise. Delay-increasing attacks can allow an adversary to distort its perceived location [6]. Delaydecreasing was also studied, e.g., by manipulating ICMP ping and traceroute as they fail to preserve the integrity of timing measurements. Combining both attacks, an adversary can forge the calculated location to an accuracy of a few tens of kilometres relative to a target desired location [7].

III. CLIENT PRESENCE VERIFICATION (CPV)
CPV [8] is a measurement-based technique designed to verify the geographic locations of web users (clients henceforth) over the Internet. The client is assumed to be motivated to misrepresent its location, e.g., to gain location-dependent benefits. CPV's design takes into consideration various adversarial location-forging tactics, including delay manipulations and IP-hiding technologies like VPNs and anonymizers. CPV does not rely fundamentally on the clients' IP addresses, nor does it determine geographic locations. Rather, it verifies an asserted (unverified) location, typically made by a client. The client's location could be asserted using, e.g., the client's GPS coordinates, the client's IP address, or even explicitly such asking the user to fill-in their street address in an online form during login.
To verify location assertions, CPV relies on an infrastructure of geographically scattered nodes, dubbed verifiers (e.g., virtual private, or cloud-based servers). The technique works as follows. When a client visits a website and asserts the geographic location from which he/she is currently browsing, three verifiers surrounding the asserted location are selected. The verifiers measure (in realtime) network delays between themselves and the client's browser, and analyze these delays to corroborate that the client is present somewhere inside the triangle determined by their (the verifiers) geographic locations. Because the verifiers cannot pinpoint where exactly the client is within the triangle, the size of the triangle is the verification granularity.

A. Secure One-way Delay Estimation
The verifiers do not measure round-trip times (RTT) between themselves and the client. They rather estimate the smaller of the forward and reverse one-way delays (OWDs) between each of them and the client. The larger OWD is discarded because propagation delays between two network nodes are bounded by the physical distance between them, so a smaller OWD measurement is a better representation to the geographic distance between both nodes than the largerthe larger must have been affected by other factors such as network congestion or circuitous routing.
To measure the OWD between a verifier and the client, CPV does not rely on standard OWD-estimation protocols like OWAMP (RFC 4656), as those require honest client cooperation, e.g., client clock synchronization and honest reporting of delays. As such, CPV relies on the Minimum-Pairs (MP) protocol [9]. MP requires the three verifiers, A, B, and C, to first synchronize their clocks, and pre-share cryptographic keys to ensure operational integrity.
Through Javascript, the client's browser is first directed to establish a WebSocket (RFC 6455) connection to the three Verifiers, which are chosen based on the client's asserted location. Verifier A begins by a sending a cryptographicallyprotected timestamp (in millisecond precision) to the client, which the browser forwards to the other two verifiers. On receiving this, verifier B calculates the propagation time from A − → client − → B, and likewise when the timestamp is received by C. Verifiers B and C then follow suit, taking turns in sending timestamps. When all three verifiers are done exchanging timestamp messages, they will have six delay values as follows: Between each pair of verifiers, e.g., between {A − → client − → B} and {B − → client − → A}, the verifiers exclude the larger OWD, and solve a system of three equations simultaneously for an estimate to the smaller OWD between the client and each verifier. That is, if the smaller of the forward and reverse OWD between the client and A, B, and C respectively is a, b, c, then (note: = sign here is used to indicate mathematical equality rather than an assignment operator): , and so on. Analysis of MP's accuracy showed that the protocol is likely to provide more accurate estimates to the smaller OWD than simply using half the RTT [9].

B. Corroborating Presence Inside the Triangle
In order to avoid potential inaccuracies from delay-todistance mapping, the calculated OWDs are not mapped to distances. They are rather compared to the smaller OWDs between the verifiers themselves, which are measured and updated periodically in a background process, independent of whether or not a client is currently being verified. Assuming x = min(AB, BA) is the smaller of the forward and reverse OWDs between verifiers A and B directly (not to be confused with min(AtB, BtA) from the previous section), and likewise y = min(BC, CB) and z = min(AC, CA), then the client's asserted location is accepted as inside the triangle if: area( xab)+area( ybc)+area( zca) ≤ area( xyz) + Such that area( xab) is the area of that triangle calculated from its side lengths x, a, and b. The value of is used to account for the two extra access network traversals occurring at the client when the timestamps propagate from a verifier to the client to another verifier (see Section III-C for setting this value).
Iterative Delay Measurement. To account for abrupt delay spikes or network irregularities, the above process of OWD calculations and comparison with those between the verifiers is iteratively repeated, e.g., 20 times. If the condition is met for the majority of the conducted iterations, the location assertion is accepted.

C. CPV Calibration
There are several parameters that tune CPV's reaction to events. The most important three are (in Section III-B), n, which is the number of delay measurement iterations, and τ , which is the fraction of those iterations that must pass (i.e., the condition is met) for the client's asserted location to be accepted. This calibration should take place before the location verification process begins. To do that, the three verifiers may use network nodes that they know as a ground truth are inside the triangle. From the network delays of these nodes, the verifiers compute values for the above mentioned three parameters, and then run CPV to verify the client's location.

D. Hindering Illicit Traffic Relaying
In an attempt to defeat geolocation, a middlebox (like a proxy server or a VPN gateway) that is physically inside the triangle can be specifically customized to filter out the verifiers' timestamps from the client's traffic and forward them to the verifiers on behalf of the client. This threat against CPV is exacerbated by the presence of numerous cheap public VPN providers whose primary service is to enable subscribers to evade geolocation technologies.
Techniques like CPV can mitigate against this by adapting known Proof-of-Work techniques [10]. The verifiers generate a cryptographic client puzzle with each timestamp message, which the client's browser must solve before forwarding the message (puzzle solution and timestamp) to the other two verifiers. The puzzles must be easy to solve so that they do not (1) overwhelm the client with high processing costs, and (2) overshadow the network propagation delays.
In the case of a middlebox connected to many simultaneous (i.e., cheating) clients, the middlebox will choose to either solve these puzzles on behalf of the clients, or forward them to the clients. In the latter case, the network delay between the middlebox and the client will get added to the time the verifiers observe for location verification, which results in CPV correctly detecting the client's absence from the respective triangle. It is thus in the middlebox's interest to choose the former case-solving the puzzles on behalf of the clients. However, this means that as more clients are connected, the middlebox will have to solve more puzzles. When these puzzles begin to accumulate, they will increase queueing delays, which contribute to the delays observed by the verifiers, eventually causing CPV to reject the location assertions of all middlebox-connected clients.
In this model, there are two main parameters contributing to the puzzle queueing rate at the middlebox: the puzzle difficulty and the middlebox's computational resources (e.g., processor clock frequency and number of processing cores). Queueing analysis [10] shows that the puzzle difficulty has a higher impact on the rate of puzzle queueing than the middlebox's computational power. This analysis suggests that this puzzle mechanism will be effective to hinder illicit middlebox relaying.

E. Evaluation Results
CPV was evaluated using PlanetLab-a distributed testbed for Internet measurement research and network experimentsusing 80 planetlab nodes in North America. Three of the nodes were selected to act as verifiers, and the remaining 77 acted as clients. Some of the 77 nodes were inside the triangle (recall: the verification triangle determined by the three verifiers), while others were outside. All 77 nodes carried out the protocol with the verifiers simultaneously to get their locations verified. Knowing the ground-truth of which nodes are inside and which are outside (the geographic locations of PlanetLab nodes are publicly disclosed on PlanetLab's website), we can count the number False Rejects (nodes inside the triangle identified by CPV as outside) and False Accepts. The process is repeated after choosing a different triangle (a different set of three nodes to act as verifiers), again counting false rejects and false accepts. In total, 34 triangles where chosen. Triangles were chosen to be nearly equilateral (physically), with inside angles ranging from 50-70 degrees (0.87-1.22 radians). The smallest triangle had an area equivalent to a circle of radius 100km, and largest 400km.
When the inside nodes are not too close to the triangle's sides (i.e., away from the closest side by at least 10% of its length), CPV resulted in a total of 1.0% false accepts and 2.0% false rejects [8]. These results were obtained when 600 CPV iterations (see Section III-B) were performed with each client. The results were not much different when only 100 iterations were performed, where the false accept rate increased only to 1.1% and the false reject rate remained unchanged. However, when only 10 iterations were performed, false accepts and false rejects were at 2.1% and 4.1% respectively.
Testing was later repeated to assess the effect of WiFi access networks on CPV's efficacy [11]. WiFi access networks often have higher delays and delay jitters. A different evaluation technique was used, as the PlanetLab infrastructure used above involved nodes connected using wired access networks. To model WiFi clients, 802.11 delay models from the literature were used to generate the last-mile delays, which were added (i.e., arithmetic summation) to the delay traces collected from PlanetLab. Since higher network delays for nodes inside the triangles may result in higher false rejects, the generated 802.11 delays were only added to the the delays of inside nodes to create the most stressful testing situation. 802.11 networks employ slotted retransmissions. The delays were generated such that each slot is 20 µsec, the propagation delay from the device to the wireless gateway is 1 µsec, and that four other wireless devices were continuously competing for the wireless media along with each wireless CPV client. With these parameters, CPV's false accepts was at 2% and false rejects at 4%. Although CPV's efficacy was affected by the WiFi access network, increasing the number of iterations can improve the results; see [11].

F. Live Demo
A live demo of CPV is currently running on http://cpv. ccsl.carleton.ca. This link hits a webserver in Ottawa, Canada, which enables clients to verify if they are present inside a USbased triangle determined by verifiers in San Francisco, Las Vegas, and San Diego. The verifiers are provided by hosting services DigitalOcean, ServerPoint, and M5 Hosting. Each VM has a 500MB RAM, and runs Ubuntu 16.04. NTP is used to synchronize their clocks. Additionally, each verifier issues an NTP query every 30 minutes using the ntpq utility to calculate the clock offset with the other two verifiers, which is added to the calculated OWDs between the verifiers for more accurate OWD estimates. Each verifier issues a timestamp to the other two verifiers every six seconds for direct OWD measurements between the verifiers.
A Java implementation of a CPV verifier runs on top of a lightweight custom-written WebSocket server, which is also implemented in Java. When a location verification request is initiated, the verifiers first check that it was issued from the authentic server (the one based in Ottawa in that demo implementation), because this server digitally signs connection IDs when they are issued. Additionally, each exchanged timestamp message between the verifiers through the client is corroborated using an MD5-based HMAC (a stronger HMAC is recommended to be used in practice). For the currently running demo, eight delay-measuring iterations are performed, once every 300ms. When all iterations are performed, the verifiers send the measured delays back to the Ottawa server, which processes the result and returns it to the browser as a jQuery response. No client puzzles are implemented yet in this demo as of this writing, nor any automatic calibration of CPV's parameters (Section III-C). Instead, the main server has manually set parameters of = 10ms and τ = 0.7, which are static and used across all clients.

IV. SERVER LOCATION VERIFICATION (SLV)
Analogous to CPV but on the server side, SLV [12] works by finding evidence of a server's physical presence inside a geographic region by measuring the server's network delays. A browser typically communicates with an SLV Manager, which orchestrates a network of server location verifiers. The challenges faced in doing so are quite different from verifying clients: (1) clients do not normally have the ability to write and run code on the server, whereas that was easily achievable by the server on the client typically using Javascript; (2) The common physical distribution of webcontent using Content Distribution Networks (CDNs) and replication technologies begs the questions: Of the multiple physical servers that may serve client content, which such servers should be selected to geograpihcally locate (verify), in order to provide a useful server-authentication service? And how should that machine be identified?
The answers to these questions depend on the threat model and the application for which geolocation is to be used. Since the goal of SLV is to reinforce server authentication, the implementation of SLV takes the view that the first machine that terminates the client's TCP (and TLS) handshake is the most critical one. The protection provided from verifying that first machine would be comparable to that provided by TLS, e.g., in the cases where the browser fetches content from multiple machines some of which are not TLS-protected (i.e., a page with mixed content).
For deciding on the mechanism used to identify machines, it is important to dissect Man-in-the-Middle (MitM) and server impersonation attacks. In MitM attacks, an adversary hijacks network traffic intended to the authentic server, and relays it to the authentic server with or without modification. Hijacking could occur on several layers of the network stack, as follows. Note that using uncompromised TLS protects against the following hijacking cases; the value of using server location to reinforce server authentication is more profound for non TLSenabled websites or to catch attacks against the TLS system.
• Case 1: Attacker's machine has a different IP address than the authentic server. In upper layers, phishing and pharming attacks are prominent traffic hijacking examples; the outbound traffic from the client has a different IP address than that of the authentic server. If the browser submits the domain name of the visited website to the SLV Manager, the Manager may resolve it to a different IP address than that seen by the browser (which could also occur benignly in the cases of CDNs). Verifying the geographic location of that IP address then becomes useless to the browser because a MitM adversary would go undetected. It is thus important to have the browser resolve a domain and submit the IP address to the SLV Manager. • Case 2: Attacker's machine has same IP address as authentic server. In lower layer hijacking (such as MAC table poisoning and ARP spoofing) and in BGP spoofing, outbound traffic from the client has the same destination IP address as that of the authentic server. Such tactics are based on routing manipulation, so that traffic intended to the authentic server's IP address reaches a different "network location" (versus geographic location), which corresponds to the attacker's machine.
In comparison to upper layers, lower layer hijacking attacks tend to be more scalable, affecting larger proportion of clients. For MAC table poisoning and ARP spoofing, the closer the attacker's machine to the authentic server's network, the more the affected clients. Likewise, BGP spoofing can cause traffic hijacking at a global scale [13]. This implies that as with higher layer traffic hijacking attacks discussed above, identifying the server by its IP address will likely allow the SLV Manager to detect if the browser-intended machine is at a different geographic location than that asserted, e.g., from a static location mapping previously obtained for that IP address.
Revisiting the above questions, this means that if SLV targets the IP address (as resolved by the browser) of the first machine that the client initially handshakes (regardless of whether the browser will be instructed to fetch other content from different places later in the session), it can detect most of the above server impersonation attacks.

A. Verification Mechanism
After obtaining an unverified server location assertion, three verifiers surrounding that location are selected. The verifiers measure network RTTs to the server over several layers, including Application using HTTP request-response times and Transport using TCP handshake responses. By means of comparing these delays with the delays between the verifiers, each pair of verifiers then verify whether the server is physically present inside the circle whose diameter is the physical distance between the verifier pair, and center is the midpoint between them (see Fig. 3).

B. Evaluation Results
Pilot testing of ∼200 experiments was conducted on SLV using PlanetLab, half of which were true location assertions made by servers and the other half were false assertions. As with CPV, the rates of false rejects and false accepts were the fundamental evaluation parameters. SLV resulted in 0% false accepts and 2.4% false rejects [12]. Although the false reject rate may seem high for some applications, it can be improved by proper selection of verifiers, e.g., those with sufficient network bandwidth and processing resources.

C. SLV Browser Extension
We have built a Firefox browser extension to reinforce TLS by integrating the webserver's verified physical location as described above into the server authentication process. The extension sends the IP address of the server to the SLV Manager, and receives the location verification result. The extension uses FlagFox to obtain an assertion (unverified) for the server's location. It also displays a flag in the URL (recall: Fig. 1), adding to that a green tick mark or a red cross indicating whether the location asserted by FlagFox is true (according to SLV's verification) or not. This process takes a few seconds to execute, during which a throbber is displayed by the flag instead. Note that such visual cues are only meant as visual feedback (in prototypes), and not an indication that we would expect end-users to take decisions upon. See below for how policies could be implemented to automatically take decisions on behalf of users.
Server location pinning in the browser. To avoid having the user interpret these icons, the SLV extension is supported with a location pinning feature, whereby a browser saves that a website (identified by its URL) was previously verified to host content from this geographic location (analogous to key pinning [14]). Although location verification is preformed based on the IP address (recall: the SLV Manager only receives an IP address from the browser), location pinning in the browser is based on the domain name. Upon receiving the verification result for a website, its location gets pinned only if the result was positive, i.e., the location was verified as true. This operation follows a trust on first use (TOFU) concept.
In general, for interpreting a received verification result, the SLV extension checks if that location (to some degree of geographic precision) was pinned before for that website. The result of the operation is threefold: Critical, Suspicious, or Unsuspicious. The first is when the verification result for a previously pinned location was negative. A Suspicious outcome is when location verification fails, but no location was pinned before for that website. Finally, an Unsuspicious outcome is when location verification passes for a domain that was not previously pinned. Note that these are only meant to illustrate how a client might utilize SLV, but we expect different applications would make different choices.
Such outcomes could result in the browser automatically taking decisions, e.g., through a policy-based engine. An example would be to instruct the browser to block/terminate the connection for all Critical outcomes of the user's personal banking website. Such a language is subject to more research scrutiny, and is not yet part of the above SLV extension.

V. CONCLUDING REMARKS
This article gave a technical overview of recent advancements in the field of secure geolocation over the Internet. Two technologies were explained, CPV and SLV, to address client and server geolocation respectively. Both rely on network timing measurements for secure location verification, taking into consideration safety measures to limit adversarial manipulations. Of the wide variety of applications that may benefit from secure location information of clients and servers, reinforcing authentication (i.e., location-aware authentication) for both ends remains an important example. Future research on CPV and SLV includes further enhancing their accuracy (in terms of the false reject and accept rates) and efficiency for large scale deployments in practice.


Tocopherols, tocotrienols and tocomonoenols: Many similar molecules but only one vitamin E

The aim of this article is to correct a very general error in scientific articles, in textbooks and in the Internet that has become an accepted fact. In this literature, the term “vitamin E″ is used for several similar molecules (both tocopherols and tocotrienols) that have never been shown to have vitamin property, i.e. a protective effect against the human deficiency disease. In fact, the name “vitamin E″ should only be used to define molecules that prevent the human deficiency disease “Ataxia with Vitamin E Deficiency” (AVED). Only one such molecule is known, α-tocopherol. This error may confuse consumers as well as medical doctors, who prescribe vitamin E without realizing that the current use of the name includes molecules of unknown, if not unwanted functions.

In addition, more recently, very similar molecules have been described, the tocomonoenols, with only one double bond in the sidechain (Fig. 2) Do all those molecules, tocopherols, tocotrienols and tocomonoenols act in the human body in a way that they can be called vitamins?
Vitamins are organic substances essential in small quantities to normal metabolism, found in minute amounts in natural foodstuffs or sometimes produced synthetically: deficiencies of vitamins produce specific disorders which can be prevented by the administration of the missing vitamin. https://www.dictionary.com/browse/vitamin. According to this definition, molecules such as essential fatty acids or essential amino acids cannot be called vitamins being needed not in minute amounts. Also trace metals, like zinc, copper, iron and molibdenum cannot be called vitamins, given their nature of inorganig compounds. In addition, vitamins are compounds whose deficiency results in a disease.
Are the molecular structures of vitamins specific? Scurvy is a disease caused by deficiency of vitamin C, historically famous because it affected sailors until the end of the 18th century, whose diets were lacking fresh fruits and vegetables. The disease can be prevented by L-ascorbic acid but not by the enantiomer D-ascorbic acid and only L-ascorbic acid can be called vitamin C (Fig. 3).
L-Ascorbic acid but not the enantiomer D-ascorbic acid prevents the deficiency disease . Only L-Ascorbic acid is called vitamin C.
Another example is given by rickets prevention by vitamin D (Fig. 4),Are all tocopherols and tocotrienols vitamins?
Tocopherols and tocotrienols act as free radical scavengers in membranes and lipoproteins [3] . They quench fatty acid peroxyl radicals by becoming tocopheroxyl radicals. Tocopheroxyl radicals potential damage would be avoided via reduction by an appropriate reducing agent (Fig. 5).
Different tocopherols have gene regulatory functions, not related to their radical scavenging capacity.
In adddition to their radical scavenging action, non-antioxidant functions of tocopherols have been shown, both in vivo and in vitro [8,9]. In T cells from mouse spleen, α-tocopherol affects signal transduction, transcriptional regulation, apoptosis pathways and other genes associated with the regulation of cell cycle. The modulation of gene expression is dependent on the fine molecular structure of α-tocopherol: S and R isomers of the molecule side chain act differently and marker genes respond uniquely to natural vitamin E (RRR-α-tocopherol) and not to the SSS isomers, that have identical radical scavenging properties [10]. The two vitamers α and γ tocopherols differ only by a methyl group on the chromane ring and have essentially the same antioxidant properties certain genes uniquely respond to either αtocopherol or γ tocopherol [11].
Different tocopherols and tocotrienols have largely different cellular effects.
β-Tocopherol has similar antioxidant properties as α-tocopherol but does not inhibit protein kinase C activity, cell proliferation and gene expression, and it is 10-fold less active than α-tocopherol in inhibiting thrombin-induced PKC activity. β-Tocopherol, has no effect on IL-1β release in monocytes, in endothelial cells, β-tocopherol, contrary to α- tocopherol has no effect on caspase-3 activity. β-Tocopherol and γ-tocopherol but not α-tocopherol strongly inhibit intracellular tyrosinase [12].
γ−Tocopherol, and its metabolite γ−CEHC, exhibit activities different from α-tocopherol such as natriuretic, anti-inflammatory, antitumoral, and more. α-Tocopherol does not exhibit anticancer properties, but it reduces anticancer actions of γ-tocopherol in vivo. δ-Tocopherol is responsible for a proinflammatory response promoted by reactive oxygen species, it prevents hormone-dependent breast cancer progression, colon carcinogenesis, lung tumorigenesis and prostate cancer cell growth. δ-Tocopherol reduces lipid accumulation, has antiangiogenic effects and elevates L-type calcium channel activity to increase neuronal differentiation [12]. In summary, the activities of β-, γ-, and δ-tocopherols do not correspond to their radical scavenging behavior: they reflect anti-inflammatory, antineoplastic, and natriuretic actions probably mediated through specific binding interactions and have no property suggesting a protective effect agains AVED.
In vitro studies with tocotrienols indicate effects against cancer (mostly due to suppression of angiogenesis), cardiovascular and neurodegenerative diseases (ability to suppress glutamate-induced activation of c-Src kinase) [2], It has been concluded that tocotrienols are more potent than tocopherols and that the different tocotrienols, which differ in their number of methyl groups and their position on the chroman ring, also differ in their biological activities. While various studies have indicated that α-tocotrienol is neuroprotective, δand γtocotrienol have been shown to exhibit the greatest anticancer effects. However, none of these molecules has been tested for an effect on AVED and only α-tocopherol has been shown to protect against AVED, the Fig. 1. The differences in the α, β, γ and δ forms are in the number and position of the methyl groups on the chromanol ring. Only the βand γforms of tocopherols or tocotrienols can be called isomers, having the same formula but a different arrangement of atoms in the molecule. The difference between the tocopherols and the tocotrienols is in the presence of three double bonds in the side chain of the latter.  [4] discovered α-tocomonoenol in palm oil, (A) and Yamamoto et al. [5] isolated an isomeric and chemically distinct α-tocomonoenol from the lipophilic fraction of salmon eggs having an unusual methylene unsaturation at the isoprenoid-chain terminus (B). only human disease associated with vitamin E deficiency.
Since it is known that small structural changes can eliminate the vitamin effect of a molecule and since only α-tocopherol, and none of the other tocopherols, have been shown to protect against the human deficiency disease, only α-tocopherol can be called bona fide vitamin E. Why the other seven fat-soluble tocopherols and tocotrienols, that have never shown any protective effect against a human deficiency disease, should be called (and sold as) vitamin E? The recently described αtocomonoenol should also be called "vitamin E"?

Conclusions
• α-Tocopherol only has been tested and shown to prevent AVED, the vitamin E deficiency disease In a more general perspective, it is not a good policy that molecules with no vitamin effect are called and sold as vitamins, especially when their unknown properties may even endanger, rather than protect, human health. For the time being, medical doctors should be aware of the confusion: if they want to prescribe vitamin E it should be clear to them that all molecules called vitamin E, except for α-tocopherol, do not have the vitamin property implied in their name. Fig. 4. Rickets is a skeletal disease caused by lack of vitamin D. Vitamin D is also structurally specific: 1,25-dihydroxycholecalciferol prevents the development of rickets but 25-hydroxycholecalciferol, just lacking one OH group (circled in the active molecule), is inactive. Consequently, only 1,25-dihydroxycholecalciferol is called Vitamin D [6].  [3] has been assumed to be at the basis of disease prevention, but this assumption was never proven. An additional claim of vitamin function for β-, γ-, and δ-tocopherol was based on their protection against fetus resorption in rats [6], but was incorrectly extended to the prevention of the human disease (AVED), the human neurological disease [7]. Based on a wrong assumption and an illogical extension of the animal effects to human disease, all tocopherols (wrongly defined frequently as "vitamin E isomers") and tocotrienols, have been therefore called vitamin E.


Improvement of the Cardiac Oscillator Based Model for the Simulation of Bundle Branch Blocks

: In this paper, we propose an improvement of the cardiac conduction system based on three modiﬁed Van der Pol oscillators. Each oscillator represents one of the components of the heart conduction system: Sino-Atrial node ( SA ), Atrio-Ventricular node ( AV ) and His–Purkinje system ( HP ). However, while SA and AV nodes can be modelled through a single oscillator, the modelling of HP by using a single oscillator is a rough simpliﬁcation of the cardiac behaviour. In fact, the HP bundle is composed of Right ( RB ) and Left Bundle ( LB ) branches that serve, respectively, the right and left ventricles. In order to describe the behaviour of each bundle branch, we build a phenomenological model based on four oscillators: SA , AV , RB and LB . For the characterization of the atrial and ventricular muscles, we used the modiﬁed FitzHugh–Nagumo ( FHN ) equations. The numerical simulation of the model has been implemented in Simulink. The simulation results show that the new model is able to reproduce the heart dynamics generating, besides the physiological signal, also the pathological rhythm in case of Right Bundle Branch Block (RBBB) and Left Bundle Branch Block (LBBB). In particular, our model is able to describe the communication interruption of the conduction system, when one of the HP bundle branches is damaged. improve the mathematical model of the heart by using four modiﬁed Van der Pol oscillators to describe the electrical activity and we modify the FitzHugh–Nagumo equations to reproduce the behaviour of the cardiac muscle cells. Each oscillator represents one of the main natural pacemakers: Sino-Atrial node Atrio-Ventricular Bundle Branch

Introduction
The human heart is a complex electro-pump that, through cycles of depolarization and repolarization, makes the propagation of the action potential, and consequently the contraction of the cardiac muscle tissue, possible. The electrical activity of the heart can be indirectly measured by the electrocardiogram (ECG), a non-invasive method based on 12 electrodes positioned on the body surface [1,2].
Considering the crucial role that the ECG signal plays in the clinical practice, different techniques were presented in the literature in order to model the dynamics of the heartbeat. The electrical and muscular activities of the heart were analyzed through both mathematical modeling and time series analysis [3]. By using mathematical modeling, it is possible to understand the electrical activation of the heart obtaining, by simulation, both normal and abnormal rhythms [4][5][6]. The generation of synthetic ECG signals with a wide range of waveform shapes and heart rates allows the modeling of the characteristics of each subpart composing the cardiac conduction system, and the understanding of its behaviour [7][8][9][10][11]. The nonlinearity and nonstationarity of the cardiovascular system make the use of nonlinear techniques for the modelling of heart activity useful [12][13][14]. Van der Pol (VdP) and Van der Mark described and modeled the behavior of heart using nonlinear relaxation oscillators [15].
FitzHugh [16] proposed an extended version of the Van der Pol equations for the generation of the action potential. In his work, he proposed a model for the emulation of the signal observed in an excitable cell of a living organism. The model is composed of two coupled, nonlinear ordinary differential equations. The first equation describes the evolution of the neuronal membrane voltage (fast action), while the second one the recovery action of the sodium and potassium channels (slow action) [17]. Grudzinski and Zebrowski modified the Van der Pol model to allow the reproduction of the time series of the action potential generated by a natural pacemaker [18,19]. In their model, the intrinsic frequency of the two pacemakers can be changed without the need to change the length of the refractory period. Although the model allows the manipulation of both the diastolic and the refractory period, the frequency obtained is too low if compared with the physiological values. Gois and Savi [20] proposed a mathematical model always based on the Van der Pol equations, but composed of three modified oscillators. This model is able to emulate the cardiac conduction system composed of a Sino-Atrial node (SA), Atrio-Ventricular node (AV) and His-Purkinje system (HS). The oscillators are connected to each other through time-delayed couplings. This model allows for generating the electrical response of the main cardiac pacemakers and obtaining the ECG wave as a composition of these signals. Although this model reproduces normal and pathological ECG signals, it does not model the activity of atrial and ventricular muscles. The characterization of muscle electrical responses proposed by Ryzhii [21] uses a quiescent excitable FitzHugh-Nagumo-type (FHN) oscillator. Successively, in [22,23], the authors improved this model by including the depolarization and repolarization waves of the atrial and ventricules generated by modified FHN systems for each ECG wave. The advanced model was able to generate normal ECG signals and several well-known rhythm disorders. In particular, it reproduces sinus tachycardia, sinus bradycardia, complete SA-AV block and complete AV-HP block [24].
In the models proposed in [20][21][22], the cardiac conduction system was treated as a network of three modified Van der Pol oscillators representing SA, AV nodes and the HP system. While SA and AV nodes can be modelled by using oscillators, the modeling of the HP system by using a single oscillator is a simplification of the cardiac activity as shown in [20][21][22] where the electrical behaviour of right and left bundle branches that composes the His-Purkinje system (Figure 1a) is modelled by using a single oscillator. By using this approach, it is not possible to characterize bundle branches of electrical conduction diseases. In the physiological condition, when the ventricular contraction occurs, the Purkinje fibers bring the impulse from the left and right bundle branches to the ventricles' myocardium. In this way, the contraction of ventricles muscle tissue and the ejection of blood outside the heart are allowed. However, the degradation of one of the bundle branches causes a defect called Bundle Branch Block (BBB) that implies an alteration of the pathways for the ventricular depolarization [25]. This work proposes a mathematical improvement of the Rizhii model [22] to better describe the heart rhythm by considering four coupled modified Van der Pol oscillators. Besides the use of two oscillators to describe SA and AV nodes, we use two additional oscillators to describe the right and left bundle branches of the His-Purkinje system. In this way, it is possible to reproduce different pathological conditions, such as the bundle branch blocks. Moreover, the FitzHugh-Nagumo (FHN) oscillators are used to describe the electrical process in the cardiac muscles. We underline that the proposed improved model corresponds to a phenomenological model [26]. The paper is organized as follows: in Section 2, the heart conduction system and mathematical models are shown. In Section 3, results of simulations are shown, and, finally, in Section 5, conclusions are drawn.

Heart Conduction System and Mathematical Models
Since the dynamic evolution of the heart action potential are close to the dynamical response of the Van der Pol oscillator (VdP) [15], this approach has been widely used in the modelling of the heart conduction system [18,20,21,[27][28][29][30][31]. In our model, we modify the heterogeneous oscillator model used for the generation of the ECG signals as proposed by Ryzhii [21]. In addition to coupled VdP oscillators for the SA and AV nodes modeling, we use two additional oscillators to model the right and left bundle branches of the HP system. Furthermore, we modify the excitable FitzHugo-Nagumo equations to take into account the two additional oscillators (right branch (RB) and left branch (LB)). Figure 2 shows the Ryzhii model and the proposed model. In order to better understand the improved model, in the first part of the following section, we provide a theoretical overview of the heart conduction system and a description of the Ryzhii model.

Conduction System of the Heart
The cardiac conduction system shown in Figure 1a is composed of a group of muscle cells characterized by their own electrical activity that generates the contraction of the cardiac muscles; the Sino-Atrial node (SA), the Atrio-Ventricular node (AV), branches of the right and left beam, and Purkinje fibers.
The SA node is placed in the right atrium and it is the natural cardiac pacemaker. These cells are known as pacemakers because they manifest a spontaneous depolarization. They are able to reach the action potential threshold faster than every other cell. Thanks to this feature, they manage the heart rhythm [32]. The depolarization rises from the SA node, spreads in the whole right and left atria and finally reaches the AV node.
Then, it proceeds along the bundle of his spreading to both left and right bundle branches depolarizing the left and right ventricles. The electrocardiogram (ECG) is a measure of how electrical activity changes during each cardiac cycle [33]. Figure 1b shows a typical ECG waveform.
The atrial depolarization phase is represented by the P wave, the ventricular depolarization by the QRS complex, and the ventricular repolarization by the T wave. In physiological conditions, the SA node works at 60-100 bpm, faster than any other cell; the only conducting path from atria to ventricles is provided by the AV node. If the SA node has a disease and it fails to generate impulses, the cells of the AV node can stimulate the heart at a rate of about 40-60 bpm [32]. If this node also presents damaged cells, the atrial-ventricular conduction block occurs. Depending on the entity of the damage, the block can be complete or partial. In this case, not all the electrical signals generated by the SA node are transmitted to the HP system. If no excitation signal is delivered from SA or AV nodes, the HP cells can fire at a rate of about 20-40 bpm. However, if the right or left bundle branch of the HP system is corrupted, the respective ventricle is not contracted.
This condition is known as Bundle Branch Block (BBB) and it implies an alteration of the pathways for ventricular depolarization. When this condition occurs, the electrical impulse can move through muscle fibers both slowing the electrical activity and changing the propagation direction of the pulses [34].

Mathematical Model
Ryzhii in [22] describes the three natural pacemakers SA, AV and the HP bundle by a system of modified VdP equations: where x i (t) and y i (t) correspond to the action potential and the transmembrane currents of the heart, are the harmonic force terms, a i > 0, u ij represent the nonlinear damping force parameters, f i are the parameters related to the intrinsic frequency of the oscillator, the coupling coefficients K SA−AV and K AV−HP represent the unidirectional coupling between the SA, AV and HP pacemakers, y τ n i ≡ y i (t − τ n ) are the velocity coupling components of the time-delay signal, and τ n are the time delays [18,20]. The synchronism between the three oscillators depends on the coupling coefficients K [35,36]. A lot of coupling methods can be found in the literature [37][38][39]. However, it is possible to avoid the use of delays τ n in the coupling terms by choosing the appropriate coupling coefficients K SA−AV and K AV−HP . In conclusion, the VdP Equations (1)-(3) become Starting from this model, we propose improvements to describe the phenomenological behaviour of the right and left bundle branch. Leaving the SA and AV node equations (Equations (4) and (5)) unaltered, we describe the right bundle branch (RB) and the left bundle branch (LB) as In the physiological condition, the right and left bundle branches are synchronized and the sum of their behaviour has to be the same as the HP oscillator: where α 1 and α 2 are coefficients used to adjust the oscillation amplitudes. In order to obtain the physiological signal, we set α 1 = 0.5 and α 2 = 0.5. In this manner, each oscillator is weighted at 50% in the generation of the signal. This is because, in a normal rhythm, the contribution of the left and the right bundle branch is the same. The parameters a i , u ij , f i , d i , e i , are chosen to obtain intrinsic oscillation rates of 70 bpm, 50 bpm, 35 bpm and 35 bpm for uncoupled SA, AV, RB and LB, respectively ( Figure 3) and with behaviour close to action potentials of real pacemakers: [40][41][42].
In particular, we use the following experimental parameters: a 1 = 40, a 2 = 50, a 3 = 50, a 1 = 40, We describe the depolarization and repolarization process for atrial (AT) and ventricular (VN) muscles starting from the modified FitzHugh-Nagumo model [16,43] proposed by Ryzhii in [22]: where k i are the scaling coefficients. In our work, we used the experimental parameters: k 1 = 2 * 10 3 , k 2 = 4 * 10 4 , k 3 = 10 4 , k 4 = 2 * 10 3 ,  The activation currents I i represent the coupling between the SA and the AT muscles and between the RB-LB pacemakers and the VN muscles. With respect to the Ryzhii model [22], we adjust the activation currents of the QRS complex and T wave in order to have the HP oscillator composed of the RB and LB oscillators: In particular, where y 3 HP = (y 3 RB + y 3 LB ), K AT De = 9 * 10 −5 , K AT Re = 1.5 * 10 −5 , K V N De = 9 * 10 −5 , K V N Re = 9 * 10 −5 , α 1 = 0.5, α 2 = 0.5 and α 3 = 0.001. In Equation (20), the terms y 3 RB and y 3 LB represent the impulses traveling through the right and left ventricles and through the bundle branches. The terms y τ RB 3 RB and y τ LB 3 LB represent the impulses traveling through the right and left ventricles through the myocardium cells. These secondary pathways are slower than the main path propagated through the bundle branches. In the physiological case, the contribution of these impulses is negligible with respect to the contribution of the main pathways [12]. On the other hand, in the pathological case, where one branch is damaged, the secondary pathways guarantee the contraction of the ventricle not served by the main pathway.
To provide a complete model, we need to take into account the contribution of the secondary pathways. It can be described as an impulse that travels through the bundle branches delayed and attenuated by a multiplicative constant α 3 . We note that, in the actual case, these secondary impulses travel through the myocardium cells and the propagation of electrical activity should be described with a specific mathematical model. However, the scope of our work is to provide a phenomenological model. For this reason, we assumed that the impulses that travel through the myocardium cells can be modeled as an impulse attenuated and delayed in time.
In the model of the left bundle branch block (LBBB), the activation currents I V N De and I V N Re represent respectively the coupling between the RB and LB pacemakers with the VN muscles; consequently, Equations (18) and (19) become For LBBB simulation, we set the experimental parameters: K SA−AV = 100, K AV−BR = 285, In the model of right bundle branch block (RBBB) the activation currents, I V N De and I V N Re are For RBBB simulation, we set the experimental parameters: K SA−AV = 100, K AV−BR = 0, K AV−LB = 285, K AT De = 9 * 10 −5 , K AT Re = 2 * 10 −5 , K V N De = 5.5 * 10 −5 , K V N Re = 7 * 10 −5 , τ RB = 0.05, α 1 = 0, α 2 = 1 and α 3 = 1.
Finally, we obtain the synthetic ECG signal as a composition of the AT and VN waveforms where z 0 = 0.2 provides the adjustment of the baseline and k R is a multiplicative coefficient to modulate the amplitude of the R peak.

Results
In order to validate our model, we performed different Simulink simulations reproducing either normal and pathological signals. Simulations are performed using a fixed step Runge-Kutta solver and a fixed step size of 0.0001 seconds. Figure 4 shows the normal condition where all pacemakers are dominated by the SA node. As a consequence, the heart rate follows the SA rate (70 bpm). In this case, both the right and left bundle branches are working, and they give the same contribution. The experimental parameters are the same used in the non-pathological as described in Section 2. Then, we performed a simulation of a pathological condition of the His bundle. We obtain this condition by changing the physiological parameters with the pathological parameters in order to obtain both complete right bundle branch block (RBBB) and complete left bundle branch block (LBBB). The pathological parameters used are shown in Section 2. Figure 5 shows the left bundle branch block (LBBB). In this case, the left ventricle is not directly activated by the impulses travelling through the left bundle branch. The right ventricle is normally activated by the right bundle branch signal. These impulses are able to travel through the myocardium of the right ventricle to the left ventricle and depolarize it (yellow arrows in Figure 5). This activation extends the QRS duration to >120 ms [32,33]. The direction of the slow depolarization (from the right to the left) produces tall R-waves in the lateral leads (I, aVL, V5-6) of the electrocardiogram and deep S-waves in the right precordial leads (V1-3). Since the ventricles are activated sequentially (first right then left) rather than simultaneously, this produces a broad or notched ('M'-shaped) R-wave in the lateral leads (I, aVL, V5-6). Figure 6a shows simulation results for the LBBB case. Experimental parameters: K SA−AV = 100, K AV−BR = 285, K AV−LB = 0, K AT De = 9 * 10 −5 , K AT Re = 3 * 10 −5 , K V N De = 9 * 10 −5 , K V N Re = −8 * 10 −5 , τ RB = 0.78, α 1 = 1, α 3 = 1 and α 2 = 0.

Left Bundle Branch Block
In this case, the term y 3 LB in Equation (21) is not present because the contribution of the left bundle branch is zero. On the other hand, the impulse generated by the right bundle brunch y 3 RB is able to both normally depolarize the right ventricle and to travel towards the left ventricle through the myocardium cells. As this second conduction way is slower than the main bundle of His-Purkinje fibres, the term y 3 RB in Equation (21) (that represents the left ventricle depolarization) is delayed through the coefficient τ RB . This produces a broad or notched ('M'-shaped) R-wave in the ECG signal (Figure 6b). It is possible to compare simulation results with a real patient's ECG signal affected by LBBB shown in Figure 6c.

Right Bundle Branch Block
In the Right Bundle Branch Block (RBBB), the activation of the right ventricle is delayed because the depolarization has to spread from the left ventricle (yellow arrows in Figure 7). The left ventricle is activated normally and the early part of the QRS complex is unchanged. The delayed right ventricular activation produces a secondary R-wave (R') in the right precordial leads (V1-3) and a wide, slurred S-wave in the lateral leads [32,33] (Figure 7).
In this case, the term y 3 RB in Equation (24) is not present because the contribution of the right bundle branch is zero. In this case, the impulse is generated by the left bundle brunch y 3 LB that is able to both normally depolarize the left ventricle and to travel towards the right one through the myocardium cells. Figure 8a shows the simulation results for the RBBB case. For this simulation, we set the following experimental parameters: K SA−AV = 100, K AV−BR = 0, K AV−LB = 285, K AT De = 9 * 10 −5 , K AT Re = 2 * 10 −5 , K V N De = 5.5 * 10 −5 , K V N Re = 7 * 10 −5 , τ LB = 0.05, α 1 = 0, α 2 = 1 and α 3 = 1. In addition, in this case, it is possible to compare simulation results with a real patient's ECG signal affected by RBBB shown in Figure 8c.

Discussion
The model discussed in this paper is an extension of the one proposed in [22]. Both the models (the original proposed by Ryzhii and the one proposed in this paper) are able to generate realistic ECG signals in order to simulate normal and pathological rhythms. The analysis of the synthetic ECG signal generated by the two models can be performed like a real signal of the patient, which is conducted by observing the morphology of different waves and their behavior in time, with specific interest on the QRS complex. This means that, for RBBB and LBBB, as for the other pathologies discussed in this paper and in [22], the anomalies are identified analyzing the shape of the signal and identifying the interest points in which the waves change their morphology. For example, if we want to reproduce the LBBB, we need to modify the shape of the QRS complex, as well as the polarity of the T wave [46]. As shown in Figure 6, selecting proper parameters, our model is able to reproduce the main typical characteristics of the signal related to this pathology.
However, it is important to note that the synthetic signal may not be exactly the same as the patient's actual ECG signal. To generate a synthetic signal closer, but not equal, to the real one, it is possible to modify, inside a specific operating range, the model parameters to obtain different nuances of the same pathological behaviour. This difference is due to the uniqueness of the signal generated by the cardiac activity of each person. It means that it is possible to compare two signals belonging to two different people by observing the waveforms in time and searching the features of interest, but not comparing the shape of the tracks point by point. In addition to personal traits, there are several factors that make any real signal different from any other-for example, the noise due to the acquisition method, the constitution of patient that influenced the transmission of the electrical signal from heart to body surface, the coexistence of different pathologies, and so on. For this reason, nowadays, the interpretation of an ECG signal is mainly performed by physicians by observing the specific characteristics of tracks.

Conclusions
In the paper, we improve the mathematical model of the heart by using four modified Van der Pol oscillators to describe the electrical activity and we modify the FitzHugh-Nagumo equations to reproduce the behaviour of the cardiac muscle cells. Each oscillator represents one of the main natural pacemakers: Sino-Atrial node (SA), Atrio-Ventricular node (AV), Right Bundle Branch (RBB) and Left Bundle Branch (LBB).
Numerical simulations show that the proposed model, based on four oscillators, is able to reproduce the heart behavior by generating the physiological signal like the model proposed by Ryzhii in [22] (based on three oscillators). The main improvement of the proposed model is the description of the His-Purkinje system with two additional oscillators, in order to better describe the behavior of each bundle branch. The model is able to reproduce the communication interruption in the heart electrical conduction system when one of the His bundle branches does not work. In this case, the model reproduces both normal ventricle contraction due to healthy bundle branch and the propagation of action potential through myocardial cells to the ventricle where the bundle branch is pathologically altered.
The waveforms obtained by simulation are comparable with the real ECG signals of patients affected by the LBBB and RBBB pathologies. In other words, the clinical phenotypes observed in the patient's signals are similar to waveform features obtained with our model. In conclusion, the proposed model makes the study of interactions between the main components of the heart possible. It allows simulation and evaluation of heart activity and dynamics under different types of pacemaker coupling. These aspects are widely useful in clinical practice.