BOUNDARY FEEDBACK STABILIZATION OF THE MONODOMAIN EQUATIONS

. Boundary feedback control for a coupled nonlinear PDE-ODE system (in the two and three dimensional cases) is studied. Particular focus is put on the monodomain equations arising in the context of cardiac electrophysiol-ogy. Neumann as well as Dirichlet based boundary control laws are obtained by an algebraic operator Riccati equation associated with the linearized system. Local exponential stability of the nonlinear closed loop system is shown by a ﬁxed-point argument. Numerical examples are given for a ﬁnite element discretization of the two dimensional monodomain equations.


1.
Introduction. This paper is concerned with the problem of boundary feedback stabilization for a coupled nonlinear reaction diffusion system. Both Neumann and Dirichlet boundary control is studied. The dynamics of interest are described by the so-called monodomain equations, a reasonably accurate simplification of the bidomain equations (see, e.g., [23]) which are frequently used to model the electrophysiological activity of the human heart. It is well-known that anomalous behavior such as fibrillation processes can be recognized within the mathematical model in the form of spiral or reentry waves. In this context, the clinical goal consists of a successful termination of these waves by external stimuli, usually called defibrillation. In [9], we have shown that the monodomain equations can locally be stabilized around a stationary solution by means of a state feedback law obtained from the linearized system. In contrast to the setup considered therein, a more realistic scenario assumes that the stimulating electrodes are situated at the boundary of the domain. With this in mind, our goal is to extend the concepts from [9] to the case of a boundary control. We further employ slightly different Lipschitz continuity estimates than the ones used in [9]. Proceeding this way, we are able to generalize the local stabilization results for the nonlinear system to a broader class of initial perturbations. Mathematically the monodomain equations are an evolution system consisting of a diffusion equation with a polynomial, more specifically, cubic nonlinearity coupled with an ordinary differential equation. Such systems play an important role far beyond their use for the description of the electrical activity of the heart. In fact, many models of this kind arise in cellular biology as described in [19], for example. One such application area is the modeling of wave propagation in excitable systems, as for instance nerve axons. They also arise in modeling chemical reaction dynamics and are known there under the name of "Schlögl" model, see e.g. [12,28].
Thus, let us consider the following nonlinear reaction diffusion system v(x, 0) =v + y 0 and w(x, 0) =w + z 0 in Ω, where Ω ⊂ R n , n ∈ {2, 3} denotes a bounded open set with smooth boundary Γ = ∂Ω and constants ι, κ > 0. The factor θ ∈ {0, 1} is used to model both the Neumann and the Dirichlet boundary control setup, respectively. The nonlinearity I ion (v, w) is assumed to be of FitzHugh-Nagumo type, i.e.
I ion (v, w) = av 3 − bv 2 + cv + dw, with a, b, c, d ∈ R + . Here, f v , f w and g are external forcing functions and (v,w) ∈ H 3 (Ω) × L ∞ (Ω), s > 0 denotes a stationary solution to 0 = ∆v − I ion (v,w) + f v in Ω 0 = ιv − κw + f w in Ω θ ∂v ∂ν = (θ − 1)v + g on Γ. ( The function m is assumed to localize the control in a part of the boundary Γ. For the precise assumptions on m, we follow the setup from [29]. In particular, if Γ is of class C 3 , we assume that m ∈ C 2 (Γ), m ≥ 0 and m(x) = 1, x ∈ Γ c , where Γ c denotes an open control subset in Γ. Similar to [9], we want to find a control input function u ∈ L 2 (0, ∞; L 2 (Γ)) such that the solution (v, w) to (1) locally decays exponentially, i.e., provided the initial perturbation y 0 = (y 0 , z 0 ) is small in an appropriate sense. For 0 ≤ σ < κ, let us introduce y = (y, z) := (e σt (v −v), e σt (w −w)), u = e σtũ and focus on where α(x) = −3av 2 + 2bv − c + σ. Instead of the stabilization of (1) around a solution of (3), we rather consider the stabilization of (4) around zero. The strategy discussed in the following is based on a linear quadratic infinite horizon optimal control problem for the linearized system and has been successfully employed for many different systems, see, e.g., [3,4,6,25,29]. In particular, for the linearized system, we closely follow the abstract Riccati theory for boundary control systems as presented in [2,21]. The structure of the paper is as follows. In Section 2 we study the case of a Neumann boundary control together with a partial Dirichlet boundary observation. Based on the general results from [21, Chapter 2] and [2], we analyze the effect of a Riccati-based feedback approach obtained for the linearized system when used in the full nonlinear setting. In particular, we slightly extend the Lipschitz continuity estimates from [9]. This will lead to a local stabilization result (in the two and three dimensional cases) for initial perturbations y 0 ∈ H ε (Ω) × L 2 (Ω), ε ∈ 1 2 , 1 . In Section 3 we investigate the possibility of using a Dirichlet boundary control. Similar to [25,29], for ε ∈ 0, 1 2 we first have to show that the involved Riccati operator Π satisfies Π ∈ L(H ε (Ω) × L 2 (Ω), H ε+2 (Ω) × L 2 (Ω)). As a result, in the two dimensional case, for ε ∈ 0, 1 2 , we obtain a local stabilization result for the nonlinear system provided the initial perturbation y 0 ∈ H ε (Ω) × L 2 (Ω) is small enough. In Section 4 we illustrate some of the theoretical results in a numerical setup. We conclude with a short summary of our contributions in Section 5.
Notation. For p ≥ 1 and s ≥ 0, we denote by L p (Ω) and H s (Ω) the usual Lebesgue and Sobolev spaces. We define H s ν (Ω) := y ∈ C ∞ (Ω) ∂y ∂ν = 0 on Γ , where C ∞ (Ω) denotes the space of infinitely differentiable functions in Ω and where the closure is taken with respect to the y H s (Ω) , s ≥ 0, norm. For s > 3 2 , we have H s ν (Ω) = y ∈ H s (Ω) ∂y ∂ν = 0 on Γ and for s ∈ [0, 3 2 ), we have H s ν (Ω) = H s (Ω), see, e.g., [8,. Furthermore, we set H s 0 (Ω) = {y ∈ C ∞ 0 (Ω)}, where the closure is taken with respect to the y H s (Ω) , s ≥ 0, norm. For s > 1 2 , we have that H s 0 (Ω) = {y ∈ H s (Ω) |y = 0 on Γ } . Given p ≥ 1, an interval I ⊂ R and a Hilbert space X, we denote with L p (I; X) (Bochner) p-integrable functions on I with values in X. We define the space For Q ∞ = Ω × (0, ∞) and r ≥ 0, s ≥ 0 we define the anisotropic Sobolev spaces H r,s (Q ∞ ) by Similar, we use Σ ∞ . We further use L 2 (Ω) instead of L 2 (Ω) × L 2 (Ω) and define In more general, boldface letters are always associated with the PDE-ODE system while italic letters refer to either the PDE or the ODE part. Let X and Y denote Hilbert spaces. For a closed densely defined linear operator A : D(A) ⊂ Y → Y, the resolvent set of A is denoted by ρ(A). Given its infinitesimal generator A with e At we denote the associated semigroup. For the definition and calculus involving interpolation spaces [X, Y] θ , we refer to e.g. [8,22,30]. By C, C 1 and C 2 we denote generic constants that may vary throughout consecutive calculations.
Hence, let us consider y(x, 0) = y 0 and z(x, 0) = z 0 in Ω, as a boundary control system as described in, e.g., [8,21,31]. For this purpose, we introduce the operators (A, D(A)) and (A * , D(A * )) with For the rest of this section, let λ 0 ∈ ρ(A) be such that A : for all y ∈ D(A), for all y ∈ D(A * ).
As a consequence, − A generates an analytic exponentially stable semigroup on Y = L 2 (Ω) and, additionally, its fractional powers A γ (see [24]) are well-defined. Further, we introduce the Neumann map N A by N A u = y iff λ 0 y − ∆y − α(x)y + dz = 0 in Ω, ∂y ∂ν = u on Γ, By explicitly solving the second equation in (7) for z, we equivalently obtain Recall that by assumptionv ∈ H 3 (Ω) which for n = 2, 3 implies thatv 2 ∈ H 3 (Ω), see, e.g., [ In particular, for s = 0 let us emphasize that , ε > 0. (9) We now may rewrite the linearized system as a boundary control system of the form d dt where B = AN A ∈ L(L 2 (Γ), [D(A * )] ) and M is the multiplication operator associated with m.
Proof. Assume that p ∈ D(A * ) and u ∈ L 2 (Γ). With y = (y, z) let us denote N A u = y. We then have . By (7) and the fact that p ∈ D(A * ), we thus obtain Finally, note that (9) implies that 2.1. Existing Riccati theory for the linearized system. With system (5), we associate the following cost functional where the multiplication operator M denotes the observable part of y on the boundary Γ and Cy = Cy where C : H 1 2 +s (Ω) → L 2 (Γ), s > 0, is the Dirichlet trace operator Cy = y |Γ . Based on the results in [2,21], our goal now is to find u in feedback form such that the cost functional (13) is minimized. We therefore make the following hypotheses which we also comment on subsequently: (H1) A is the infinitesimal generator of a strongly continuous analytic semigroup on Y := L 2 (Ω).
(H5) There exists an operator K ∈ L(Z, D( A δ )), such that the strongly continuous analytic semigroup e (A+K M C)t on D( A δ ), generated by A + K M C is exponentially stable.

TOBIAS BREITEN AND KARL KUNISCH
As we have shown in [9, Lemma 3.1], A indeed is the infinitesimal generator of a strongly continuous analytic semigroup on Y. Hence, (H1) is satisfied. From (12), it follows that (H2) holds for U = L 2 (Γ) and Lemma 2.1 implies that C = B * and M C ∈ L(D( A δ ), Z) for such that (H3) is satisfied. From now on, let us fix 1 4 < δ < 3 8 . It remains to show the validity of (H4) and (H5). For this, we may follow the strategy used in [9] and first consider the decoupled problem where ω > 0, together with the cost functional which is finite forỹ ∈ L 2 (Q ∞ ) and u ∈ L 2 (Σ ∞ ). Using available null controllability results from [15, Theorem 2] and a well-known extension argument (see, e.g., [16, Theorem 2.3]) we particularly conclude that for each y 0 ∈ L 2 (Ω) there exists u ∈ L 2 (Σ ∞ ) such that the corresponding solutionỹ satisfiesỹ ∈ L 2 (Q ∞ ) so that J(u,ỹ) < ∞. In other words, the decoupled system (16) is stabilizable. As a consequence, there exists a linear feedback operator K ∈ L(L 2 (Ω), L 2 (Γ)) such that the decoupled closed-loop operator A := ∆ + α + ω + BK generates an exponentially stable semigroup on L 2 (Ω). According to [9, Lemma 3.2], if we choose ω = ε + ιd (κ−σ)−ε , 0 < ε < κ − σ, by means of a Schur-complement type argument, one can show that all eigenvalues of the operator are located in the open left complex half plane. This implies that the coupled system (5) is in fact ε-stabilizable and assumption (H4) is fulfilled. Finally, the detectability condition (H5) of the pair (A, M C) is equivalent to the stabilizability of the pair (A * , C * M ). Since we already know that B = C * , the arguments provided for (H4) apply here as well and imply (H5). In summary, all of the assumptions (H1)-(H5) are indeed fulfilled for the monodomain equations and we can use the results from [21, Section 2.5] that guarantee the existence of a unique nonnegative self-adjoint solution Π = Π * ∈ L(Y) to the algebraic Riccati equation for all z 1 , z 2 ∈ D( A δ+s ), s > 0 : Additionally, by [21, Theorem 2.5.2], it holds that and for the transformed system whereȳ = A δ y and the modified cost functional In particular, for the corresponding Riccati operatorΠ it is shown that and that the semigroup generated by 2.2. The nonhomogeneous linear system. Before we turn to the nonlinear setting, we first consider the nonhomogeneous equation with is an isomorphism.
The idea now is to apply this isomorphism to the (transformed) system (23) associated with AΠ. In particular, for AΠ, we have to characterize the interpolation spaces [D(AΠ), Y] α and [D(A * Π ), Y] 1−α . As pointed out in [2], the fact that A has bounded imaginary powers together with the perturbation results from [14, Proposition 2.7] allows us to conclude that AΠ has bounded imaginary powers.
According to [30,Theorem 1.15 Following [2], let us then introduce H s Π := D(A s Π ) and recall from [2] the characterization of H s Π : In summary, Theorem 2.2 implies that is an isomorphism. Note that for a result like (25), in [2, Corollary 2] the author assumes thatC satisfiesC * C ∈ L(D(( A * ) 2 )) which does not hold in our case. Before we proceed with system (23), we give a more precise characterization of H s Π , see also [2, Section 3/4] for the corresponding results in the context of the Navier-Stokes equations.
As a consequence of the previous results, for (23), we finally obtain.
Proof. Due to the relation between Π andΠ, it holds that Hence, instead of (23), let us consider the transformed system 2 )] which, due to (24), again means that ). Once more, with Proposition 1, we find that Since z = A δ y, from Proposition 1 we find that  1 2 ).

2.3.
Results for the nonlinear system. In [9] we have shown some Lipschitz estimates for nonlinearities of the form (2). In what follows, we extend these results to cover a broader range of parameters ε. .
Proposition 2 allows to show the following results.
For the nonlinearity in (4), i.e., we obtain the following local Lipschitz continuity property.
In particular, we have that .

TOBIAS BREITEN AND KARL KUNISCH
Proof. Following the arguments provided for similar results in [9,26], let us show that the mapping M : z → y z defined by is a contraction in D µ . As in [9], define with C 1 and C 2 as in Theorem 2.3 and Lemma 2.5, respectively.
Step 1. For z ∈ D µ , from Lemma 2.5 we get Utilizing Theorem 2.3 we obtain that Hence, M is mapping D µ to itself.
Step 2. For z 1 = z As a consequence, from Theorem 2.3 we know that Using Lemma 2.5, this yields The mapping M is a contraction in D µ and the proof is complete.
As a consequence, for system (1), we obtain the following result.
With λ 0 and A as in Section 2 let us define the Dirichlet map D A via D A u = y iff λ 0 y − ∆y − α(x)y + dz = 0 in Ω, y = u on Γ, As before, by explicitly solving the second equation for z, and using well-known results [ Consequently, for the linearized system we now obtain the following boundary control system d dt where B = AD A ∈ L(L 2 (Γ), [D(A * )] ) and M is the multiplication operator associated with m.
Proof. For p ∈ D(A * ), u ∈ L 2 (Γ) and D A u = y we have .
By (29) and the fact that p ∈ D(A * ), we thus obtain Similar as before, let us note that 3.1. Riccati theory for the linearized system. Let us consider the following cost functional corresponding to (27) Obviously, for the observation operator Cy(t) = y(t) it holds that C ∈ L(Y, Z), where Z = L 2 (Ω). As in the Neumann case, conditions (H1)-(H5) hold. From [21, Section 2.1] we now get the existence of a unique nonnegative self-adjoint solution Π = Π * to an algebraic Riccati equation. We further know ([21, Theorem 2.2.1]) that ( A * ) ν Π ∈ L(Y) for any 0 ≤ ν < 1. As emphasized in [21] and investigated in [5,25], in certain cases it is even possible to take ν = 1 such that Π ∈ L(Y, D( A * )).
In particular, in [2] it was shown that this always holds true in case that C * C ∈ L(D( A 2 )). Note that the latter condition is fulfilled for C as in (34). Moreover, for the case of the Navier-Stokes equations, in [25], it has been further shown that Π ∈ L H ε (Ω), H 2+ε (Ω) ∩ H 1 0 (Ω) , ε ∈ 0, 1 2 . With regard to the desired local stabilization result for the nonlinear system, a property similar to the latter one will be essential.

3.2.
Results for the nonlinear system. As in the case of Neumann boundary conditions, we first consider a nonhomogeneous equation of the form with f = f 1 0 .
With Theorem 2.2 and Proposition 3, we obtain the following result.
Let Ω ⊂ R 2 and ε ∈ 0, 1 2 . Then there exist µ 0 > 0 and a nondecreasing function η from R + into itself such that if µ ∈ (0, µ 0 ) and y 0 H ε (Ω)×L 2 (Ω) ≤ η(µ), then admits a unique solution in the set Remark 1. For the results in Theorem 2.6 and Theorem 3.3 we assumed that U = L 2 (Γ). However, in practical applications it is often more appropriate to utilize finite dimensional controllers. In this case, one typically assumes that the control function u is separable, i.e., u( where h i are shape functions and denotes the dimension of the control space U = R . An essential tool for showing stabilizability by finite dimensional controllers then is the decomposition of the spectrum into an infinite dimensional stable part and a finite dimensional unstable part, see, e.g., [13,Chapter 5]. In essence, the stabilizability problem then reduces to its finite dimensional counterpart such that, given shape functions h i , i = 1, . . . , , one may utilize the Hautus test for stabilizability. Let us emphasize that due to the structure of the monodomain equations, in [10], we have explicitly specified the spectrum in terms of the spectrum of the PDE part. In particular, as long as the desired stabilization rate is limited by the stability of the ODE, the unstable part of the system is discrete. This way, one may apply the exact same arguments to extend the results in Theorem 2.6 and Theorem 3.3 to the case of finite dimensional controllers. We end this section with our final result for the local exponential stabilization of system (1) in case that θ = 0.

Numerical examples.
For the numerical validation of the theory presented in this manuscript, we study the following version of the monodomain equations where Ω = (0, 1) × (0, 1) and all other parameters are to be specified below. A finite element discretization is obtained by the MATLAB R PDE toolbox. All results correspond to a 64 × 64 regular grid with n = 2 · 4225 = 8450 degrees of freedom. Given a stationary solution of the form (3), the finite dimensional approximation of (4) with σ = 0 is given by where A n , E n ∈ R n×n and B ∈ R n× . The nonlinearity F n : R n → R n is defined elementwise according to the nonlinearity in (37). For the discrete control operator, denotes the dimension of the finite dimensional control space (cf. Remark 1). By C n ∈ R p×n we denote the finite dimensional approximation of the output operator appearing in (13) and (34), respectively. In both Neumann and Dirichlet case, the multiplication operator m is taken from [29]. The closed loop system is obtained by solving the following generalized algebraic matrix Riccati equation For this, a Kleinman-Newton iteration as described in, e.g., [7,11,20], has been used. As a stabilizing initial guess Π (a) The Neumann case.  In our experiments, we chose (v 3 ,w 3 ) for which the corresponding linearized system exhibited an unstable subspace of dimension 2.
Perturbation around constant stationary state. In Figure 2, numerical results for an initial value of the form (v 0 , w 0 ) = (v 3 + ξ v ,w 3 + ξ w ) with ξ v = 0.3 · randn(n,1) and ξ w = 0.003 · randn(n,1) are given. As is seen in Figure 2a the uncontrolled system remains close to the unstable stationary for a certain period of time before it slowly starts to decrease in magnitude. Finally, for t = 1000, the system approaches the stable stationary solution (v 1 ,w 1 ) from below. On the other hand, Figure  2b demonstrates the successful local stabilization of the nonlinear system. The different dynamical behavior between uncontrolled and controlled system is even more evident from the results shown in Figure 2d. Here, the temporal evolution of the L 2 (Ω) error between computed and desired state is visualized. Again we emphasize that the uncontrolled system remains close to the desired state first. Figure 2c shows the control law for some of the control domains specified in Figure  1a. In particular, we see that the magnitude of the feedback control is relatively small when compared to the actual state.
Stabilization of a reentry wave. As a second test case, we study the feedback law when the initial value (v 0 , w 0 ) is chosen such that it causes a reentry pattern (as arising in context of fibrillation processes). To be more precise, in Figure 3a we   illustrate the wave-like dynamics appearing for the uncontrolled system. Snapshots of the evolution for the closed loop system are provided in Figure 3b. We emphasize that, as we have also seen in the first example, when the control is switched off close to the stationary solution (v 3 ,w 3 ), the system will converge to the stable resting state (v 1 ,w 1 ). This way, the feedback law implicitly allows to perform a defibrillation process by first controlling the system to the unstable stationary state. From there the uncontrolled dynamics converges to the origin. Again, a more detailed comparison between controlled and uncontrolled dynamics can be obtained in terms of the L 2 (Ω) error, see Figure 3d. For a precise pattern of the control laws, we refer to 2c. In order to ensure that (v 3 ,w 3 ) is also a stationary solution of the PDE system, in (37) we assume an inhomogeneous Dirichlet boundary condition g =v 3 . The resulting linearized system has 6 eigenvalues with positive real part. In contrast to the Neumann case, we only take three controls that are visualized in Figure 1b.   Perturbation around constant stationary state. The results shown in Figure 4 correspond to a similar setup as in the Neumann case. For both, uncontrolled and controlled system, the initial value was taken as a perturbation (0.3 randn(n,1)) of the unstable stationary solution (v 3 ,w 3 ). As a consequence, the uncontrolled system shows an oscillatory behavior, see Figure 4a. On the other hand, the feedback law computed from the solution of the algebraic Riccati equation allows to locally stabilize the nonlinear system, see Figure 4b. For a more detailed comparison, we refer to the temporal evolution of the L 2 (Ω) error that is provided in Figure 4d. The performance of the three boundary control patches is visualized in Figure 4c, also underlining the convergence of the closed loop system.

5.
Conclusions. We have studied boundary feedback control problems for a class of nonlinear reaction diffusion systems of PDE-ODE type. A particular emphasis was on the so-called monodomain equations. For the linearized system, we have investigated the use of classical Riccati-based feedback controllers. Based on some new Lipschitz estimates for the cubic non monotone nonlinearity, we could show that the feedback law derived for the linearized system locally stabilizes the nonlinear system as well. While in the case of Neumann boundary conditions, both the two and three dimensional case could be handled, in the Dirichlet case we restricted