An adjoint approach to identification in electromyography: modeling and first order optimality conditions

We start with a mathematical description of the human body. To avoid computational overhead, we consider only a locally truncated part, containing the region of interest. This is a reasonable assumption, since static electric fields decay quickly far away from the source.

To represent some part of the body geometrically, we choose a bounded open domain ${\Omega}\subset {\mathbb{R}}^{3}$ with Lipschitz boundary ∂Ω. That means locally, ∂Ω is the graph of a Lipschitz continuous function. In [8, 31] a more detailed definition is given. Since the human body consists of different tissue types, we split the domain Ω into several subdomains Ω_i , such that Ω = ∪Ω_i . For simplicity's sake, we restrict our model to three different tissue types, namely muscle, fat, and bone tissue. We use Ω_M (muscle tissue), Ω_F (fat tissue), and Ω_B (bone tissue) to denote the corresponding subdomains of Ω, which are all assumed to have Lipschitz boundary, as well.

Figure 3 shows such a described domain which represents some part of the hand. Here the domain is truncated at the fingers and at the wrist to reduce the numerical effort. From the truncation, we get an artificial boundary, which we call ∂Ω_A. Finally, skin tissue bounds the rest of our domain, and we label it with ∂Ω_S. We suppose that the skin part has positive boundary measure and

$$\begin{equation*}\partial {\Omega}=\partial {{\Omega}}_{\text{S}}\cup \partial {{\Omega}}_{\text{A}},\qquad \varnothing =\partial {{\Omega}}_{\text{S}}\cap \partial {{\Omega}}_{\text{A}}.\end{equation*}$$

Each tissue type is equipped with three different electromagnetic properties: conductivity σ, permittivity , and permeability μ. Since human tissue is not magnetic its permeability is given by the permeability of vacuum μ₀. We also know that the upper frequency limit in human tissue is around 1 kHz, see [22, 23, 27]. We assume that the tissue is purely resistive in this frequency range, and thus the properties are independent of the exact frequency, see [22, 30]. Furthermore, it is well known that bone and fat tissue are isotropic, concerning conductivity and the permittivity. In contrast, muscle tissue is anisotropic, where conductivity and permittivity are higher in the direction of the muscle fibers, see [1, 22, 30]. Hence, properties depend on the geometry of the muscle. But for simple geometries, we may assume that the muscle fibers are straight or only slightly curved. Thus we can represent the conductivity tensor by a 3 × 3 matrix which is constant in each muscle.

Tabular 1 lists the conductivity and permittivity values for different tissue types. Permittivity is given relatively to the permittivity of vacuum, which is ${{\epsilon}}_{0}=8.8e-12\enspace \text{A}\enspace \text{s}\enspace {\text{V}}^{-1}\enspace {\text{m}}^{-1}$. The conductivity is positive and constant in these tissue type and thus the two estimates

$$\begin{equation}{\sigma }^{\mathrm{max}}{:=}\mathrm{max}\left\{{\Vert}{\sigma }^{M}{{\Vert}}_{\infty },\vert {\sigma }^{F}\vert ,\vert {\sigma }^{B}\vert \right\}< \infty ,\end{equation} \tag{ 2 }$$

$$\begin{equation}{\sigma }^{\mathrm{min}}{:=}\mathrm{min}\left\{{\lambda }_{\mathrm{min}}(\sigma ),\vert {\sigma }^{F}\vert ,\vert {\sigma }^{B}\vert \right\} > 0,\end{equation} \tag{ 3 }$$

are valid. This shows that the conductivity is in L^∞(Ω) and elliptic, which is needed for ellipticity and boundedness of the bilinear form in (1).

Finally, we discuss the a geometric model of a motor unit, which is a bundle of muscle fibers. Figure 1 shows the sketch of a motor unit. Motor units are the smallest entities within a muscle that can be controlled individually by the brain. First, our brain sends a signal to the neuromuscular junction, which lies in the innervation zone of the muscle. The innervation zone lies thereby in the middle of the muscle. When a motor unit is activated, it generates two action potentials that propagate toward both ends of the motor unit, see [14] and cause a contraction of the muscle cells along the way. These propagating action potentials also create an electric potential in the whole tissue, which can be measured on the skin.

Zoom In Zoom Out Reset image size

It is known, that all muscle fibers in a motor unit are activated almost simultaneously, consequently only a single fibered motor unit is measured, see [22]. On the one hand, this complicates the inverse problem since different fiber configurations of a motor unit can cause similar measurements. On the other hand, we can simplify the simulation by simulating only a single fibered motor unit. Figure 1 sketches that the neuromuscular junctions and the ends of the fibers are uniformly distributed in small areas of the motor unit, see also [14, 22]. When compared with the motor unit, the extension of such an area is negligible. Furthermore, by linearity of Poisson's equation (1), we can apply the superposition principle, i.e., we can sum up the action potentials of each muscle fiber to one action potential and scale he action potential (13) of a single representative fiber by the approximate number of muscle fibers in the motor unit. Thus in non-pathological cases, it is legitimate to use a single fibered motor unit for the inverse problem. If needed, we can then conclude from the single fibered motor unit the distribution of the muscle fibers by using statistical methods. We represent such a single fibered motor unit by a regular curve $u\in {H}^{2}(-1,1,{\mathbb{R}}^{3})$ and denote by u(0) the neuromuscular junction and by u(−1) and u(1) the two ends of the motor unit.

A surface EMG measurement device measures the electric potential at the skin that is caused by a moving electric charge within the motor unit. To simulate a measurement we need a model that connects a moving electric charge with the corresponding electric potential in the tissue. The velocity of these moving charges is relatively slow, such that electrodynamic effects (like emission of electromagnetic waves) can be neglected. Thus, we will use a quasi-static model, see [22, 27].

We denote the moving electric charge by ρ(x, t), and the corresponding potential by Φ(x, t), where x is the spatial variable and t the time. Then, in classical form we obtain the following model, which would have to be augmented by boundary conditions and transmission conditions to take into account jumps of σ at material boundaries:

$$\begin{equation}-\mathrm{div}\enspace \sigma (x)\nabla {\Phi}(x,t)=\rho (x,t),\end{equation} \tag{ 4 }$$

see [27]. Moreover, we will model the action potential as a moving charge along a fiber, which is a curve in 3D. So ρ(⋅, t) is rather a measure than a function in x. We will thus proceed to derive an appropriate weak form of (4).

As mentioned above, there are two different boundary types, namely ∂Ω_S, which represents the skin, and ∂Ω_A, which is an artificial boundary, caused by the truncation of the domain. For ∂Ω_S we assume that there are no other sources outside of the domain, and thus the potential is zero there. Those assumptions lead to the following Robin boundary condition

$$\begin{equation*}\sigma (x)\nabla {\Phi}(s,t)\cdot \nu +\mu {\Phi}(s,t)=0\quad \;\text{at}\;\enspace \partial {{\Omega}}_{\text{S}},\end{equation*}$$

where μ > 0 is the skin conductivity. For ∂Ω_A, we also assume that no other sources are present in the truncated part of the body. That means we can use homogeneous Neumann boundary conditions in this case:

$$\begin{equation*}\sigma (x)\nabla {\Phi}(s,t)\cdot \nu =0\quad \;\text{at}\;\enspace \partial {{\Omega}}_{\text{A}}.\end{equation*}$$

Combining these aspects, formal integration by parts yields the following symmetric bilinear form

$$\begin{equation}\begin{aligned}\hfill a& :{H}^{1}({\Omega})\times {H}^{1}({\Omega})\to \mathbb{R}\hfill \\ \hfill a({\Phi},v)& {:=}\underset{{\Omega}}{\int }(\sigma (x)\nabla {\Phi}(x))\cdot \nabla v(x)\enspace \mathrm{d}x+\underset{\partial {{\Omega}}_{\text{S}}}{\int }\mu {\Phi}(s)v(s)\enspace \mathrm{d}s.\hfill \end{aligned}\end{equation} \tag{ 5 }$$

Due to (2), (3), and the presence of Robin boundary conditions with μ > 0 on ∂Ω_S, this bilinear form is H¹(Ω)-elliptic by a generalized Poincare inequality, cf e.g., [31, lemma 2.5]. Thus, by the Lax–Milgram theorem, we obtain a continuously invertible linear operator

$$\begin{align*}\hfill A:{H}^{1}({\Omega})\to {H}^{1}{({\Omega})}^{\ast }\\ \hfill (A{\Phi})(v){:=}a({\Phi},v).\end{align*}$$

Since H¹(Ω) is reflexive, we may identify H¹(Ω)^** ≅ H¹(Ω) and also consider the adjoint operator A*: H¹(Ω) → H¹(Ω)* of A as (A*v)(Φ) = a(Φ, v) = (AΦ)(v).

As already mentioned, we will model the electric charge at a time-instant t by a measure $\rho (t)\in \mathcal{M}(\bar{{\Omega}})$. Here $\mathcal{M}(\bar{{\Omega}})$ is the Banach space of Radon measures on $\bar{{\Omega}}$ which is isomorphic to the dual $C{(\bar{{\Omega}})}^{\ast }$ of the space of continuous functions by the well known Riesz representation theorem. Thus, we may introduce the following weak form:

$$\begin{equation}a({\Phi}(t),v)={\int }_{\bar{{\Omega}}}v\enspace \mathrm{d}\rho (t)\quad \forall \enspace v\in {C}^{\infty }(\bar{{\Omega}}).\end{equation} \tag{ 6 }$$

Since ${H}^{1}({\Omega})\;/ \hookrightarrow C(\bar{{\Omega}})$ in our 3D-setting, we cannot write (6) as an operator equation AΦ(t) = ρ(t) in H¹(Ω), and thus the Lax–Milgram theorem cannot be applied directly. Nevertheless, by an approach due to Stampacchia, see [29], solvability of (6) with Φ(t) ∈ ${W}^{1,{p}^{\prime }}$(Ω) for some p' < 3/2 can be established.

In this approach the bilinear form (5) is redefined on different spaces as

$$\begin{equation*}{a}_{p}:{W}^{1,{p}^{\prime }}({\Omega})\times {W}^{1,p}({\Omega})\to \mathbb{R},\end{equation*}$$

with 1/p + 1/p' = 1 for p > 3, implying that ${W}^{1,p}({\Omega})\hookrightarrow C(\bar{{\Omega}})$. This gives rise to the following restricted pre-dual problem for some l ∈ ${W}^{1,{p}^{\prime }}$(Ω)* ↪ H¹(Ω)*:

$$\begin{equation*}\;\text{find}\;\psi \in {W}^{1,p}({\Omega}):{a}_{p}(v,\psi )=l(v)\quad \forall \enspace v\in {H}^{1}({\Omega}).\end{equation*}$$

By Lax–Milgram, this problem has a solution ψ ∈ H¹(Ω), and it is a question of regularity theory, if ψ is an element of W^1,p (Ω). If this is true for all l ∈ ${W}^{1,{p}^{\prime }}$(Ω)*, which is known as 'maximal regularity', then the pre-dual operator

$$\begin{align*}\hfill {}^{\ast }A_{p}:{W}^{1,p}({\Omega})& \to {W}^{1,{p}^{\prime }}{({\Omega})}^{\ast }\hfill \\ \hfill ({}^{\ast }A_{p}\psi )(v)& ={a}_{p}(v,\psi ),\hfill \end{align*}$$

is an isomorphism by the open mapping theorem.

Remark 1. More generally, it can be shown that ψ is an element of ${H}^{1}({\Omega})\cap C(\bar{{\Omega}})$ if l ∈ ${W}^{1,{p}^{\prime }}$(Ω)* for p' > d/(d − 1), ${\Omega}\subset {\mathbb{R}}^{d}$, cf e.g. [17]. Then, with some additional technical effort, we can still show solvability of (6), but an additional criterion is required to single out a unique solution. A detailed discussion can be found in [25].

For simplicity we thus impose the following assumption:

Assumption 2.1. The domain Ω and its subdomains Ω_j are sufficiently regular, such that the operator ${}^{\ast }A_{p}:{W}^{1,p}({\Omega})\to {W}^{1,{p}^{\prime }}{({\Omega})}^{\ast }$ is an isomorphism for some p > 3.

Under this assumption and by reflexivity of Sobolev spaces, we can conclude that the adjoint of ${A}_{p}{:=}{({}^{\ast }A_{p})}^{\ast }$

$$\begin{align*}\hfill {A}_{p}:{W}^{1,{p}^{\prime }}({\Omega})\to {W}^{1,p}{({\Omega})}^{\ast }\\ \hfill ({A}_{p}\phi )(w){:=}{a}_{p}(\phi ,w),\end{align*}$$

is also an isomorphism, since adjoints of isomorphisms in normed spaces are isomorphisms, as well.

Due to the continuous and dense embedding ${W}^{1,p}({\Omega})\hookrightarrow C(\bar{{\Omega}})$ we can use the corresponding adjoint embedding $C{(\bar{{\Omega}})}^{\ast }\hookrightarrow {W}^{1,p}{({\Omega})}^{\ast }$ to regard the charge ρ(t) as an element of W^1,p (Ω)* for each t, and we obtain unique solvability of the operator equation:

$$\begin{equation*}{A}_{p}{\Phi}(t)=\rho (t).\end{equation*}$$

Hence, a unique electric potential Φ(t) ∈ ${W}^{1,{p}^{\prime }}$(Ω) that satisfies (6) exists for each ρ(t). Since all spaces are reflexive, we can identify the adjoint operator and the pre-adjoint operator ${A}_{p}^{\ast }={}^{\ast }A_{p}$.

With the above model, the potential Φ(t) can in principle be computed in the whole domain for every t if ρ(t) is given. However, the computational effort to do so with finite elements is too large, given that we are only interested in a certain number of measurements y_i (t) at the boundary of Ω. We thus develop a more efficient adjoint approach to compute a desired measurement $y(t)\in \mathbb{R}$ from given ρ(t).

In our setting, the potential is measured with small circular electrodes on the skin as follows

$$\begin{equation*}y(t){:=}B({\Phi}(t))=\frac{1}{\vert D\vert }\underset{D}{\int }{\Phi}(s,t)\enspace \mathrm{d}s,\end{equation*}$$

where D ⊂ ∂Ω_S is the area of the electrode. The trace theorem, cf e.g. [31, theorem 2.1], implies that B is well defined as an element of ${W}^{1,{p}^{\prime }}$(Ω)*.

Let ρ(t) ∈ W^1,p (Ω)* and denote by Φ(t) ∈ ${W}^{1,{p}^{\prime }}$(Ω) the solution of

$$\begin{equation}({A}_{p}{\Phi}(t))(v)=\rho (t)(v)\quad \forall \enspace v\in {W}^{1,p}({\Omega}).\end{equation} \tag{ 7 }$$

Now consider the solution ω ∈ W^1,p (Ω) of the adjoint problem:

$$\begin{equation}({A}_{p}^{\ast }\omega )(\phi )=B(\phi )\quad \forall \enspace \phi \in {W}^{1,{p}^{\prime }}({\Omega}),\end{equation} \tag{ 8 }$$

which corresponds to the following Poisson problem in weak form:

$$\begin{equation}\underset{{\Omega}}{\int }(\sigma (x)\nabla \phi (x))\cdot \nabla \omega (x)\enspace \mathrm{d}x+\underset{\partial {{\Omega}}_{S}}{\int }\mu \phi (s)\omega (s)\enspace \mathrm{d}s=\frac{1}{\vert D\vert }\underset{D}{\int }\phi (s,t)\enspace \mathrm{d}s\quad \forall \enspace \phi \in {W}^{1,{p}^{\prime }}({\Omega}).\end{equation} \tag{ 9 }$$

Then we compute easily

$$\begin{equation}y(t)=B({\Phi}(t))=({A}_{p}^{\ast }\omega )({\Phi}(t))={a}_{p}({\Phi}(t),\omega )=({A}_{p}{\Phi}(t))(\omega )=\rho (t)(\omega ).\end{equation} \tag{ 10 }$$

That means we can compute the potential at an electrode efficiently by evaluating

$$\begin{equation}y(t)=\underset{{\Omega}}{\int }\omega (x)\enspace \mathrm{d}\rho (t),\end{equation} \tag{ 11 }$$

where $\omega \in {W}^{1,p}({\Omega})\hookrightarrow C(\bar{{\Omega}})$ is the previously computed solution of the adjoint problem (8). In the following section, we will give a physiologically meaningful definition of the measure ρ(t), concentrated on a curve, such that (11) can be evaluated as a line integral. Then the computation of y(t) requires just the evaluation of this line integral, which is much cheaper than computing the solution of an elliptic equation. For electrodes B_i , i ∈ 1...n_E , we obtain y_i (t) via solutions ω_i of the corresponding problems ${A}_{p}^{\ast }{\omega }_{i}={B}_{i}$.

For our forward problem (6), we may assume that the moving charge is completely contained in the muscular subdomain Ω_M, which is disjoint with the domains of measurement D_i ⊂ ∂Ω_s. We can thus invoke regularity results to obtain more smoothness of the restriction $\omega {\vert }_{{{\Omega}}_{\text{M}}}$. This is useful to render sensitivities of y(t) with respect to perturbations of the support of ρ(t) well defined, which in turn is needed for the later defined optimal control problem (19).

Lemma 2.2. The solution ω of (8) is in C^∞(Ω_M) ∩ W^1,p (Ω).

Proof. By (9) ω is the solution of an elliptic equation subject to inhomogenous Robin boundary conditions on D, but without interior source terms.

Let us consider the restriction $\omega {\vert }_{{{\Omega}}_{\text{M}}}$ to Ω_M. We observe that $\omega {\vert }_{{{\Omega}}_{\text{M}}}$ satisfies a homogenous Laplace equation in weak form on Ω_M subject to Dirichlet boundary conditions on ∂Ω_M, given simply by the condition that $\omega {\vert }_{\partial {{\Omega}}_{\text{M}}}$ is the trace of ω on ∂Ω_M. Moreover, the coefficient σ is constant on Ω_M. Such problems, however, are known to be C^∞ regular in the interior, cf e.g. [13, corollary 8.11]. □

In this section we provide a model for the charge ρ(t) that moves along a fiber inside the motor unit. As described in section 2.1, the source density ρ is given by two action potentials moving along a motor unit. For simplicity, we assume that the measured data corresponds only to one active motor unit. This is a reasonable assumption, since it is possible to identify the activity of single motor units through EMG decomposition methods, see for example [19].

Since the radius of the fiber is very small, we will model our fiber as (the trace of) a fixed curve u in Ω_M. The action potential extends spatially along that fiber, but also moves along the fiber during the time span of the activation. The action potential takes a characteristic shape, sketched in figure 2, and we have to map this signal onto a certain time varying segment of the curve that describes the fiber to obtain the line measure ρ(t) at time t. Additional difficulties arise toward the ends of the fiber. Here the principle of conservation of charges is violated and therefore additional source terms are added.

Zoom In Zoom Out Reset image size

Biomedical modeling of action potentials starts with the following function

$$\begin{equation}{i}_{m}(z){:=}\begin{cases}-c({\sigma }_{\text{in}},r,{n}_{F})\mathrm{exp}(az)\left(6az+6{(az)}^{2}+{(az)}^{3}\right)\quad \hfill & \quad \text{if}\enspace z\leqslant 0\hfill \\ 0\quad \hfill & \quad \text{else,}\hfill \end{cases}\end{equation} \tag{ 12 }$$

in terms of a reference parameter $z\in \mathbb{R}$. Here a > 0 is a scaling factor that determines the spatial extension of the signal and c(σ_in, r, n_F ) is a constant depending on the intracellular conductivity σ_in, the radius of the motor unit r, and the number of muscle fibers n_F . For a more detailed description of the action potential we refer to [1, 22, 28]. The antiderivative of i_m is

$$\begin{equation}{I}_{m}(z){:=}\begin{cases}-c({\sigma }_{\text{in}},r,{n}_{F})\mathrm{exp}(az)\left(3{(az)}^{2}+{(az)}^{3}\right)\frac{1}{a}\quad \hfill & \text{if}\enspace z\leqslant 0\hfill \\ 0\quad \hfill & \text{else},\hfill \end{cases}\end{equation} \tag{ 13 }$$

and thus

$$\begin{equation*}{\int }_{\mathbb{R}}{i}_{m}(z)\enspace \enspace \mathrm{d}z={I}_{m}(0)-\underset{z\to -\infty }{\mathrm{lim}}\enspace {I}_{m}(z)=0,\end{equation*}$$

which corresponds to the principle of conservation of charge in the body. Up to now, the action potential is defined as a function on $\mathbb{R}$, so the next step is to define a pull-pack of i_m onto the given curve u ∈ Ω_M. A common assumption in biomedical modeling is that the velocity ν with which the action potential propagates along the fiber is constant, see [22]. Since the curve u represents the trajectory of the two propagating action potentials, we therefore choose the parameterization of the curve u such that it matches with the propagation velocity ν of the signal, in other words $\vert \dot{u}(\tau )\vert \equiv \nu$. That means we can identify each point on the curve u with some $z\in \mathbb{R}$ via the arc length

$$\begin{equation*}z(\tau )=\underset{0}{\overset{\tau }{\int }}\vert \dot{u}(\xi )\vert \enspace \mathrm{d}\xi =\underset{0}{\overset{\tau }{\int }}\nu \enspace \enspace \mathrm{d}\xi =\nu \tau .\end{equation*}$$

That means z(0) = 0 corresponds to the neuromuscular junction u(0) and we can identify points that are on the 'right' side of the neuromuscular junction with some $z\in {\mathbb{R}}^{+}$ and points on the 'left' side with some $z\in {\mathbb{R}}^{-}$. To model the action potential that propagates from the neuromuscular junction toward the 'right' end of the fiber, we shift the origin of the action potential i_m (z) by ν ⋅ (t₀ − t) and set

$$\begin{equation*}\tilde{\rho }(u(\tau ),t)=\begin{cases}{i}_{m}\left(z(\tau )+\nu \cdot ({t}_{0}-t)\left. \right)\right)\quad \hfill & \quad \text{if}\enspace z(\tau ) > 0\hfill \\ 0\quad \hfill & \;\quad \text{else.}\hfill \end{cases}\end{equation*}$$

To model the second action potential that propagates in the opposite direction, we mirror the signal at point zero, which is equal to adding a minus sign before z(τ) and changing z(τ) > 0 to z(τ) < 0. By combining both action potentials we get a line-measure ρ_l as follows:

$$\begin{align}\hfill {\rho }_{\text{l}}(u(\tau ),t){:=}{\hat{\rho }}_{\text{l}}(\tau ,t)& {:=}\left\{\begin{cases}\hfill {i}_{m}\left(z(\tau )+\nu \cdot ({t}_{0}-t)\left. \right)\right)\hfill & \hfill \quad \text{if}\enspace z(\tau )< 0,\hfill \\ \hfill {i}_{m}\left(-z(\tau )+\nu \cdot ({t}_{0}-t)\left. \right)\right)\hfill & \hfill \quad \text{if}\enspace z(\tau ) > 0\hfill \end{cases}\right\}\hfill \\ \hfill & ={i}_{m}\left(\nu \cdot (\vert \tau \vert +{t}_{0}-t)\right).\hfill \end{align} \tag{ 14 }$$

Conservation of charge and end-effects. As we have observed, our model i_m of the action potential respects conservation of charge: its total integral over $\mathbb{R}$ vanishes. However, our definition of y(t) involves only an integral over a bounded subset of $\mathbb{R}$ and the corresponding total charge is given by (taking into account the substitution of variables formula):

$$\begin{align}\hfill {\rho }_{\mathrm{t}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}(t)& =\underset{-1}{\overset{0}{\int }}{i}_{m}\left(\nu \cdot (-\tau +{t}_{0}-t)\right)\nu \enspace \enspace \mathrm{d}\tau +\underset{0}{\overset{1}{\int }}{i}_{m}\left(\nu \cdot (\tau +{t}_{0}-t)\right)\nu \enspace \enspace \mathrm{d}\tau \hfill \\ \hfill & ={\left.-{I}_{m}(\nu \cdot (-\tau +{t}_{0}-t))\right\vert }_{-1}^{0}+{I}_{m}{\left.(\nu \cdot (\tau +{t}_{0}-t))\right\vert }_{0}^{1},\hfill \end{align} \tag{ 15 }$$

which is non-zero in general. This truncation of the integral, can be observed in figure 2 when comparing the prototype of the action potential (see figure 2(a)) with the two propagating action potentials in the figures 2(b) and (d). In the beginning the tail parts of the action potentials are not present on the motor unit (see figure 2(b)). After some the first part of the action potentials is not longer on the motor unit (see figure 2(d)). This truncation needs to be compensated for. Otherwise, the principle of conservation of charge would be violated and yield characteristic artifacts in simulations, so called end-effects. Those end effects can be observed in the simulation shown in red (see figure 4(c)).

The representation of (15) by boundary terms at τ = −1, 0, 1 already suggests how to construct an appropriate compensation. We will add point (Dirac) measures at u(0), u(1) and u(−1), scaled by the negatives of the corresponding boundary terms to the measure ρ(t). This can also be interpreted physiologically: at the ends of the fibers, transitional imbalances of charge are compensated by small displacements of charge in the close vicinity of the end-plates. A similar approach can be found in [14].

Figure 2(b) shows that the first truncation of the action potential is at the neuromuscular junction. In (15) this corresponds to the boundary term at τ = 0. Therefore we get

$$\begin{equation*}{\rho }_{\text{s}}(u(0),t){:=}2{I}_{m}(\nu \cdot ({t}_{0}-t)),\end{equation*}$$

which is a Dirac measure at the neuromuscular junction u(0). Figure 2(c) shows that after some time the support of the action potential lies on the motor unit, which means that the imbalance of charges coming from the truncation at the neuromuscular junction tends to zero, exponentially.

When the action potentials arrive at the ends u(−1) and u(1) of the motor unit, see action potentials in figure 2(d), they are again truncated. Thus we have to compensate the boundary terms at τ = −1 and τ = 1 in (15) by Dirac measures at u(−1) and u(1) as follows:

$$\begin{equation*}{\rho }_{\text{s}}(u(-1),t)={\rho }_{\text{s}}(u(1),t){:=}-{I}_{m}(\nu \cdot (1+{t}_{0}-t)).\end{equation*}$$

Observe that these charges are 0, if t ⩽ 1 + t₀, since I_m (z) = 0 for z ⩾ 0. Figure 4(c) shows a comparison of simulations of y(t) with and without the source term compensation. Without compensation the simulated signal has two characteristic extra spikes (shown in red). These effects can be very pronounced, compared to the remaining signal, due to the smoothing effect of the potential equation. Our above described compensation technique can eliminate these end-effects, as seen in figure 4(c) by observing the green function.

Finally, to define the simulated measurement y(u, t) for a single measuring electrode, where we stress the dependence of the measurement y on u by including u as an argument, we insert ρ_l and ρ_s into the definition of y(u, t) in (11). This yields the following compensated model:

$$\begin{align}\hfill y(u,t)& ={\int }_{{\Omega}}\omega (x)d({\rho }_{\text{l}}+{\rho }_{\text{s}})(t)\hfill \\ \hfill & =\underset{-1}{\overset{1}{\int }}\omega (u(\tau )){\rho }_{\text{l}}(u(\tau ),t)\vert \dot{u}(\tau )\vert \enspace \mathrm{d}\tau +\sum\limits _{\tau \in \left\{-1,0,1\right\}}\omega (u(\tau )){\rho }_{\text{s}}(u(\tau ),t)\hfill \\ \hfill & =\underset{-1}{\overset{1}{\int }}\omega (u(\tau ))\nu {\rho }_{\text{l}}(u(\tau ),t)\enspace \mathrm{d}\tau +\sum\limits _{\tau \in \left\{-1,0,1\right\}}\omega (u(\tau )){\rho }_{\text{s}}(u(\tau ),t)\hfill \\ \hfill & =\underset{-1}{\overset{1}{\int }}\omega (u(\tau ))\nu {\hat{\rho }}_{\text{l}}(\tau ,t)\enspace \mathrm{d}\tau +2\omega (u(0))\nu {I}_{m}(\nu \cdot ({t}_{0}-t))\hfill \\ \hfill & \quad -\left(\omega (u(-1))+\omega (u(1))\right){I}_{m}(\nu \cdot ({t}_{1}-t)).\hfill \end{align} \tag{ 16 }$$

Clearly, y(u, t) also depends on the solution ω of (8) and thus on the domain D if the corresponding electrode, used for the measurement. If, for i = 1...n_E , electrodes D_i are considered, we denote the corresponding measurements by y_i (u, t).

To illustrate the properties of the forward problem, we perform the simulation of a surface EMG measurement, using the previously established model. We consider a measurement for a single fibered motor unit in the first dorsal interosseous (FDI) muscle of the right hand. The maximal extension of the hand is approximately 15 cm from the wrist to the fingers, 10 cm from the little finger to the thumb, and between 2 and 5 cm from the back of the hand to the front of the hand. At this point, we would like to thank the authors of [26] for sharing the STL files of their MRI measurements.

For simplicity, our model contains only the FDI muscle and the first two metacarpal bones. The rest of the domain was modeled as fat tissue. Figure 3 shows the geometrical model and a grid of 24 circular electrodes (white dots) placed above the FDI muscle. The location of the motor unit is depicted by the black straight line that crosses the electrode grid horizontally. Its depth below the electrode grid is 4 mm. The motor unit is represented by a piecewise cubic Hermite polynomial on 20 subintervals, which would also allow the representation of curved motor units to high accuracy.

Zoom In Zoom Out Reset image size

To incorporate the electrode grid into the STL files, we used the CAD software Blender [4]. To generate a mesh from the STL data, we used gmsh [12]. We performed all the following computations in C++, where we used the toolbox Dune [3] for all mesh-related operations, and the finite element toolbox Kaskade7 [15] to compute the finite element discretization of the adjoint problem.

The numerical computation of the adjoint solutions is done by a finite element method on a triangulation $\mathcal{T}$ of Ω consisting of 414 195 tetrahedra. On $\mathcal{T}$ we used continuous piecewise quadratic ansatz functions to discretize W^1,p (Ω) and ${W}^{1,{p}^{\prime }}$(Ω) by

$$\begin{equation*}{W}_{h}{:=}\left\{w\in C({\Omega},\mathbb{R}):w{\vert }_{K}\in {P}_{2}(K)\quad \forall \enspace K\in \mathcal{T}\right\}.\end{equation*}$$

A Galerkin method, applied to the adjoint problem (8) leads to the discrete problem

$$\begin{align*}\hfill & \text{find}\;\;{\omega }_{h}\in {W}_{h}\;\text{s.t.}\;\hfill \\ \hfill & \qquad ({A}_{p}^{\ast }{\omega }_{h})(\phi )=B(\phi )\quad \forall \enspace \phi \in {W}_{h}.\hfill \end{align*}$$

After finite element discretization, we end up with a large sparse linear system of equations to be solved. We used the preconditioned conjugated gradient method from the linear algebra library Eigen [11] to solve this linear system. Since there is no grid hierarchy available we used an standard incomplete Cholesky decomposition, see [20], as a preconditioner.

The evaluation of the line integrals (16) is performed by numerical quadrature along the motor unit, i.e., the trajectory of u. As seen in figure 2, the action potential is only nonzero on a small part of the trajectory but shows large oscillations there. Thus, a standard piecewise quadrature rule on uniform intervals would be inefficient. We therefore use an adaptive multigrid quadrature algorithm, cf e.g., [6], using Gauss–Kronrod quadrature formulas.

To evaluate the finite element function ω_h (x) at a quadrature point x ∈ Ω, we need to know the tetrahedron $K\in \mathcal{T}$ that contains x. To find K efficiently, we exploit that the quadrature points are ordered along the trajectory, and therefore we can use a neighborhood search: if x_i ∈ K_j is known, and x_i+1 is the next quadrature point, we test all neighbors of K_i if they contain x_i+1. If that fails, we find, by successive bisection, a point $\tilde{x}$ on the line segment [x_i , x_i+1] for which a neighborhood search is successful: $\tilde{x}\in \tilde{K}$. Then we repeat the whole procedure with (x_i , K_i ) replaced by $(\tilde{x},\tilde{K})$, until a neighborhood search for x_i+1 is successful. If this algorithm fails, or if an initial inclusion x ∈ K is not known, we fall back to a hierarchic search over the whole grid.

Figure 4(a) visualizes the simulated time dependent signal on all 24 electrodes. Depending on their location the electrodes yield different measurements. For example, the strength of the measured signal depends on the distance of the electrodes and the source.

Zoom In Zoom Out Reset image size

The identification of the depth of a source from boundary measurements is often difficult. We thus perform a variation in depth of the motor unit and compare the simulation result. Figure 4(b) shows the simulated measurement of one electrode (marked in red in figure 4(a)) for motor units with different depths. The simulated measurements show that most parts of the simulated potential decrease very quickly if the depth is increased. But due to the concentrated stationary sources, the potential decrease is much slower at the end of the measurement. That is a well known effect when modeling monopolar signals, see [10, 14].

As we will discuss in the following, our adjoint approach reduces the required numerical effort per simulated measurement significantly, at least for high temporal resolution. For example, the computation of 200 time steps by a direct approach would require the solution of 200 PDEs of the form (1). By an adjoint approach, we only have to solve the PDE (8) 24 times, which is the number of electrodes used. In this example the numerical effort concerning the solution of PDEs is reduced by a factor of 8.

Once, the weighing functions are computed, the numerical cost of the quadrature formulas, needed for simulation with our adjoint approach is almost negligible. In the above described example we observed the following computational times: on a standard workstation the single threaded computation of one PDE solution by a cg iteration requires about 7 s (assembly of the problem data and setup of the preconditioner not included), while the simulation of all 24 measurements with the adjoint approach required 0.3 s for all 200 time steps in total, and thus around 1.5 milliseconds per time step.

Certainly, these results depend on the spatial and temporal resolution of the problem, and also parallelization is possible in both the direct and the adjoint approach, but the numbers give a clear impression of the advantages of our adjoint approach, already for a single simulation. If multiple simulations with the same geometry have to be performed, for example inside an optimization algorithm, the weighing functions only have to be computed once. Then the computational savings of the adjoint approach, compared to a direct simulation are even more pronounced.

To prove that the problem has at least one solution, we first need some auxiliary results. First we show that the admissible set is weakly closed and that the equality constraint satisfies enough regularity. The second result shows that the functional J is differentiable and weakly lower continuous. We use then this two results two show that the objective functional satisfies the properties that we need to proof the existence result.

Lemma 3.1. The admissible set U_ad is weakly closed and the equality constraint G is Fréchet differentiable.

Proof. To show that U_ad is weakly closed, it is sufficient to show that the set

$$\begin{equation*}{U}_{1}=\left\{v\in {H}^{2}(-1,1,{\mathbb{R}}^{3})\vert v(\tau )\in {\bar{{\Omega}}}_{\text{M}}\enspace \text{for}\;\;\text{a. e}\;\;\tau \in [-1,1]\right\},\end{equation*}$$

and

$$\begin{equation*}{U}_{2}=\left\{v\in {H}^{2}(-1,1,{\mathbb{R}}^{3})\vert G(v)(\tau )=0\enspace \text{for}\;\;\text{a.}\;\text{e.}\;\;\tau \in [-1,1]\right\},\end{equation*}$$

are weakly closed. The admissible set U_ad is then also weakly closed as an intersection of finitely many weakly closed sets.

Let $u:[-1,1]{\mapsto}{\bar{{\Omega}}}_{\text{M}}$ be a regular curve. From [5, theorem 8.8] we know that a compact embedding ${E}_{1}:{H}^{2}(-1,1,{\mathbb{R}}^{3}){\mapsto}C(-1,1,{\mathbb{R}}^{3})$ exist, and thus a $v\in {H}^{2}(-1,1,{\mathbb{R}}^{3})$ exists such that E₁ v = u. Thus U₁ is well defined and not empty.

Let now {u_k } ⊂ U₁ be a weak convergent sequence with limit u. Since the embedding E₁ is compact there exist a subsequence ${u}_{{k}_{l}}$ such that ${E}_{1}{u}_{{k}_{l}}\to {E}_{1}u$ in $C(-1,1,\mathbb{R})$. Furthermore, there exist an another subsequence ${E}_{1}{u}_{{k}_{{l}_{i}}}$ that converges pointwise to E₁ u for all τ ∈ [−1, 1]. Since $\bar{{\Omega}}$ is closed and we can conclude from ${E}_{1}{u}_{{k}_{{l}_{i}}}(\tau )\in \bar{{\Omega}}$ that ${E}_{1}u(\tau )\in \bar{{\Omega}}$ and thus also u(τ). This shows that U₁ is weakly closed.

As above there exits a compact embedding ${E}_{2}:{H}^{1}(-1,1,{\mathbb{R}}^{3}){\mapsto}{W}^{1,4}(-1,1,{\mathbb{R}}^{3})$ with ${E}_{2}\dot{u}=\dot{u}$. We can conclude from Hölder's inequality that ${\langle \dot{u},\dot{u}\rangle }_{2}\in {H}^{1}(-1,1,\mathbb{R})$ and thus

$$\begin{align*}\hfill & \hfill G:{H}^{2}(-1,1,{\mathbb{R}}^{3}){\mapsto}{L}^{2}(-1,1,\mathbb{R})\\ \hfill & \hfill [G(u)](\tau )={\langle \dot{u}(\tau ),\dot{u}(\tau )\rangle }_{2}-{\nu }^{2},\end{align*}$$

is well defined. It is well known that, as a continuous bilinear form, G is Fréchet differentiable with derivative

$$\begin{equation*}[{G}^{\prime }(u)(v)](\tau )={\langle \dot{u}(\tau ),\dot{v}(\tau )\rangle }_{2}.\end{equation*}$$

Let now {u_k } ⊂ U₂ be a weakly convergent sequence with limit u. As above there exist a subsequence ${u}_{{k}_{l}}$ such that $E{u}_{{k}_{l}}\to Eu$ in ${L}^{4}(-1,1,{\mathbb{R}}^{3})$. Since G is differentiable it is also continuous. This means that

$$\begin{equation*}0=G(E{u}_{{k}_{l}})\to G(Eu).\end{equation*}$$

□

Remark 3. The H²-regularization term is essential for weak closedness of U₂. It is not hard to construct a zig-zagging sequence of trajectories u_n with $\vert {\dot{u}}_{n}\vert =\nu$ a.e. (and thus bounded in W^1,∞), such that the weak limit $\bar{u}$ violates $\vert \dot{\bar{u}}\vert =\nu$. From a computational point of view, such a regularization term does not impose severe difficulties, since u can be easily discretized in an H² conformal fashion by a piecewise polynomial spline which is globally C¹.

Lemma 3.2. The functional $J:{H}^{2}(-1,1,{\mathbb{R}}^{3})\supset {U}_{\text{ad}}{\mapsto}\mathbb{R}$ is continuous and weakly lower semi-continuous. If u(τ) ∈ Ω_M for all τ ∈ [−1, 1], then J is Fréchet differentiable at u.

Proof. First we note that, as a sum, J is continuous, Fréchet differentiable and weakly lower semi-continuous, if J₁, J₂ and J₃ are continuous, Fréchet differentiable and weakly lower semi-continuous.

To show the three properties for J₁, we define for fixed t and k the mapping

$$\begin{align*}\hfill {\psi }_{k}& :{{\Omega}}_{\text{M}}\times \left[-1,1\right]{\mapsto}\mathbb{R}\hfill \\ \hfill {\psi }_{k}(x,\tau )& {:=}{\omega }_{k}(x)\nu {\hat{\rho }}_{\text{l}}(\tau ,t),\hfill \end{align*}$$

with derivative

$$\begin{align*}\hfill {\psi }_{k,x}& :{{\Omega}}_{\text{M}}\times \left[-1,1\right]{\mapsto}\mathcal{L}({\mathbb{R}}^{3},\mathbb{R})\hfill \\ \hfill {\psi }_{k,x}(x,\tau )v& {:=}{\langle \nabla {\omega }_{k}(x),v\rangle }_{2}\nu {\hat{\rho }}_{\text{l}}(\tau ,t).\hfill \end{align*}$$

From lemma 2.2 we now that ${\omega }_{k}:{{\Omega}}_{\text{M}}{\mapsto}\mathbb{R}$ is in C^∞(Ω_M) and therefore also ∇ω_k ∈ C^∞(Ω_M). We can conclude now that ψ_k (⋅, τ) and ψ_k,x (⋅, τ) are continuous for all τ ∈ [−1, 1]. Thus ψ_k is Fréchet differentiable, see [35, page 192]. Furthermore, there exists a continuous compact embedding $E:{H}^{2}(-1,1,{\mathbb{R}}^{3}){\mapsto}C(-1,1,{\mathbb{R}}^{3})$ with Eu = u, see [5, theorem 8.8]. Therefore, the superposition operator

$$\begin{align*}\hfill & {{\Psi}}_{k}:{U}_{\text{ad}}{\mapsto}C(-1,1,{\mathbb{R}}^{3})\hfill \\ \hfill & ({{\Psi}}_{k}(u))(\tau ){:=}{\psi }_{k}(Eu(\tau ),\tau ),\hfill \end{align*}$$

is well defined and Fréchet differentiable, see [2, theorems 6.3 and 6.7]. The Fréchet derivative is given by

$$\begin{align*}\hfill & {{\Psi}}_{k}^{\prime }:{U}_{\text{ad}}{\mapsto}\mathcal{L}\left({H}^{2}(-1,1,{\mathbb{R}}^{3});C(-1,1,\mathbb{R})\right)\hfill \\ \hfill & ({{\Psi}}_{k}^{\prime }(u)v)(\tau )={\psi }_{k,x}(Eu(\tau ),\tau )Ev(\tau ).\hfill \end{align*}$$

By the same argumentation the correction terms ω_k (u(⋅))I_m (⋅) are well defined and Fréchet differentiable. Thus,

$$\begin{align*}\hfill {y}_{k}(u,t)& =\hspace{-1pt}\underset{-1}{\overset{1}{\int }}{{\Psi}}_{k}(u)(\tau )\enspace \mathrm{d}\tau +2{\omega }_{k}(u(0)){I}_{m}(\nu \cdot ({t}_{0}-t))-\left({\omega }_{k}(u(-1))+{\omega }_{k}(u(1))\right){I}_{m}(\nu \cdot ({t}_{1}-t)),\hfill \end{align*}$$

is well defined and Fréchet differentiable with derivative

$$\begin{align*}\hfill {D}_{u}{y}_{k}(u,t)(v)& =\underset{-1}{\overset{1}{\int }}\langle \nabla {\omega }_{k}(u(\tau )),v{(\tau )\rangle }_{2}\nu {\hat{\rho }}_{\text{l}}(\tau ,t)\enspace \mathrm{d}\tau +2\langle \nabla \omega (u(0)),v{(0)\rangle }_{2}{I}_{m}(\nu \cdot ({t}_{0}-t))\hfill \\ \hfill & \quad -\left(\langle \nabla {\omega }_{k}(u(-1)),v{(-1)\rangle }_{2}+\langle \nabla {\omega }_{k}(u(1)),v{(1)\rangle }_{2}\right){I}_{m}(\nu \cdot ({t}_{1}-t)).\hfill \end{align*}$$

The chain rule indicates that ${J}_{1}(u)=\frac{1}{2}{\int }_{0}^{T}{\Vert}y(u,t)-{y}_{m}(t){{\Vert}}_{{n}_{E},2}^{2}\enspace \mathrm{d}t$ is Fréchet differentiable with derivative ${J}_{1}^{\prime }(u)(v)=\underset{0}{\overset{T}{\int }}{\langle y(u,t)-{y}_{m}(t),{D}_{u}y(u,t)(v)\rangle }_{{n}_{E},2}\enspace \mathrm{d}t$.

Finally, since ${J}_{1}:C(-1,1,{\mathbb{R}}^{3}){\mapsto}C(-1,1,{\mathbb{R}}^{3})$ is differentiable, it is also continuous and thus lower semi-continuous in $C(-1,1,{\mathbb{R}}^{3})$. By the compact embedding $E:{H}^{2}(-1,1,{\mathbb{R}}^{3}){\mapsto}C(-1,1,{\mathbb{R}}^{3})$ we conclude that J₁ is weakly lower semi-continuous in ${H}^{2}(-1,1,{\mathbb{R}}^{3})$.

J₂ and J₃ are obviously convex quadratic bilinear forms and it is well known that they are Fréchet differentiable from L₂ into itself and weakly lower semi-continuous, cf e.g. [9]. The derivatives of J₂ and J₃ are then given through

$$\begin{align*}\hfill {J}_{2}^{\prime }(u)v& =\underset{-1}{\overset{1}{\int }}{\alpha }_{1}\langle u(\tau )-{u}_{\text{ref}}(\tau ),v{(\tau )\rangle }_{2}\enspace \mathrm{d}\tau \quad \;\text{and}\hfill \\ \hfill {J}_{3}^{\prime }(u)v& =\underset{-1}{\overset{1}{\int }}{\alpha }_{2}\langle \ddot{u}(\tau ),\ddot{v}{(\tau )\rangle }_{2}\enspace \mathrm{d}\tau .\hfill \end{align*}$$

The sum rule indicates that J is continuous and Fréchet differentiable and since the lim inf is super-additive, J is weakly lower semi-continuous. □

Collecting the derivatives from the previous proof, we get for the derivative J

$$\begin{align*}\hfill {J}^{\prime }(u)v& =\underset{0}{\overset{T}{\int }}{\langle y(u,t)-{y}_{m}(t),{D}_{u}y(u,t)(v)\rangle }_{{n}_{E},2}\enspace \mathrm{d}t\hfill \\ \hfill & \quad +\underset{-1}{\overset{1}{\int }}{\alpha }_{1}{\langle u(\tau )-{u}_{\text{ref}}(\tau ),v(\tau )\rangle }_{2}+{\alpha }_{2}{\langle \ddot{u}(\tau ),\ddot{v}(\tau )\rangle }_{2}\enspace \mathrm{d}\tau .\hfill \end{align*}$$

Remark 4. Recall that the derivatives ∇ω(x), required in the definition of D_u y(u, t) are well defined by lemma 2.2, since u is assumed to be contained in the muscular tissue Ω_M. However, if ω is approximated by a finite element function, ∇ω is only piecewise continuous.

Lemma 3.3. The objective functional F is weakly lower semi-continuous and radially unbounded.

Proof. From lemma 3.2 we already now that J is weakly lower semi-continuous. It remains to show that the indicator function is lower semi-continuous.

A function $f:X{\mapsto}\mathbb{R}$ is weakly lower semi-continuous if the level sets N_α f = {x ∈ X|f(x) ⩽ α} are weakly closed for all $\alpha \in \mathbb{R}$, see [7, theorem 7.4.11]. For the indicator function the level sets are given through

$$\begin{equation*}{N}_{\alpha }{\iota }_{{U}_{\text{ad}}}(u)=\begin{cases}{U}_{\text{ad}}\hfill & \hfill \text{if}\enspace \alpha > 0\quad \\ \varnothing \hfill & \hfill \text{else.}\quad \end{cases}\end{equation*}$$

From lemma 3.1 we know that U_ad is weakly closed and since the empty set is always weakly closed we can conclude that the indicator function is weakly lower semi-continuous.

Obviously we have that F(u) > 0 and for ||u||_2,2 → ∞ either J₂, J₃ or ${\iota }_{{U}_{\text{ad}}}$ goes to infinity and therefore F is radially unbounded. □

These three auxiliary results enable us to prove existence of optimal solutions by standard techniques:

Theorem 3.4. The optimization problem (20) has at least one solution in ${H}^{2}([-1,1],{\mathbb{R}}^{3})$.

Proof. Let {u_n } be a minimizing sequence. Since F is radially unbounded by lemma 3.3 u_n is bounded. Thus by reflexivity of ${H}^{2}(-1,1,{\mathbb{R}}^{3})$ there exist a weakly convergent subsequence ${u}_{{n}_{k}}$ with limit u_*. By lemma 3.3 F is weakly lower semi-continuous and thus

$$\begin{equation*}\mathrm{inf}\enspace F\leqslant F({u}_{\ast })\leqslant \underset{k\to \infty }{\mathrm{liminf}}\enspace F({u}_{{n}_{k}})=\mathrm{inf}\enspace F,\end{equation*}$$

which shows that the limit point u_* is a minimizer of F. □

Next we derive first-order optimality conditions for the problem. This is complicated by the geometric constraint $u\subset {\bar{{\Omega}}}_{\text{M}}$, which we imposed to assert existence of solution. However, from a practical point of view, we can expect that the optimal solution u_* is contained in Ω_M without the need to enforce this as a constraint, because the measured values originate form a true signal, emitted from u ⊂ Ω_M. In addition, the part J₂ of the objective functional is penalizes by J₂ the departure of u from the reference trajectory u_ref, which reasonably will be chosen to lie in Ω_M. Therefore, for simplicity, we will drop these geometric constraints from now on and assume the optimal solutions u_* lies in the muscle domain Ω_M. That means we can rewrite the optimal control problem (19) as a purely equality constrained problem:

$$\begin{equation}\begin{aligned}\hfill & \underset{u\in {H}^{2}(-1,1,{\mathbb{R}}^{3})}{\mathrm{min}}J(u)\hfill \\ \hfill & \;\text{s.t.}\;\quad G(u)(\tau )=0\quad \;\text{for}\;\;\text{a.}\;\text{e.}\;\enspace \tau \in [-1,1].\hfill \end{aligned}\end{equation} \tag{ 21 }$$

As usual we eliminate the equality constraint with the help of a Lagrange multiplier, which lead to the following result:

Theorem 3.5. Let u_* be a local minimizer of (21) that lies in Ω_M. Then there exist a Lagrange multiplier $\lambda \in {H}^{1}(-1,1,\mathbb{R})$, such that

$$\begin{equation}\begin{aligned}\hfill 0& ={J}^{\prime }({u}_{\ast })+{G}^{\prime }{({u}_{\ast })}^{\ast }\lambda \hfill \\ \hfill 0& =G({u}_{\ast })(\tau )\quad \;\text{for}\;\;\text{a.}\;\text{e.}\;\enspace \tau \in [-1,1]\hfill \\ \hfill \lambda & \in {H}^{1}{(-1,1,\mathbb{R})}^{\ast }\hfill \end{aligned}.\end{equation} \tag{ 22 }$$

Proof. First, we recall that $\lambda \in {H}^{1}{(-1,1,\mathbb{R})}^{\ast }$ is called Lagrange multiplier, if

$$\begin{equation}\lambda \in {K}^{+},\end{equation} \tag{ 23 }$$

$$\begin{equation}\lambda \left(G({u}_{\ast })\right)=0,\end{equation} \tag{ 24 }$$

$$\begin{equation}{J}^{\prime }({u}_{\ast })-{G}^{\prime }{({u}_{\ast })}^{\ast }\lambda \in C{({u}_{\ast })}^{+},\end{equation} \tag{ 25 }$$

is fulfilled, cf [36, equation (1.1)]. Here K is a convex closed cone such that G(u) ∈ K and K⁺ is the dual cone of K. Since we have pure equality constraints, $K=\left\{0\right\}\subset {H}^{1}(-1,1,\mathbb{R})$ and from the definition of the dual cone we can conclude that ${K}^{+}={H}^{1}{(-1,1,\mathbb{R})}^{\ast }$, see [31, 36]. $\left[G({u}_{\ast })\right](\tau )=0$ implies that condition (24) is always fulfilled and thus we can replace it by $\left[G({u}_{\ast })\right](\tau )=0$. Furthermore, C(u_*) is the conical hull of ${H}^{2}(-1,1,{\mathbb{R}}^{3})$, which is ${H}^{2}(-1,1,{\mathbb{R}}^{3})$ and therefore the dual cone $C{({u}_{\ast })}^{+}=\left\{0\right\}$. Thus, condition (25) simplifies to

$$\begin{equation*}{J}^{\prime }({u}_{\ast })-{G}^{\prime }{({u}_{\ast })}^{\ast }\lambda =0.\end{equation*}$$

Since we already know from lemma 3.2 that J and G are Fréchet differentiable, we can conclude from [31, theorem 6.3] that a Lagrange multiplier λ exists, if G satisfies the regularity condition of Zowe and Kurcyusz, which is given through

$$\begin{equation*}{G}^{\prime }({u}_{\ast })C({u}_{\ast })+K(-G({u}_{\ast }))={H}^{1}(-1,1,\mathbb{R}).\end{equation*}$$

Here K(⋅) is the convex hull of K. We can conclude from K = {0} and G(u_*) = 0 that K(G(u_*)) = {0} and since $C({u}_{\ast })={H}^{2}(-1,1,{\mathbb{R}}^{3})$, this condition is equivalent to G'(u_*) surjective. To show this, we choose for arbitrary $w\in {H}^{1}(-1,1,\mathbb{R})$

$$\begin{equation*}v(\tau )=\underset{-1}{\overset{\tau }{\int }}\frac{{\dot{u}}_{\ast }(\xi )w(\xi )}{{\nu }^{2}}\enspace \enspace \mathrm{d}\xi ,\end{equation*}$$

which is in ${H}^{2}(-1,1,{\mathbb{R}}^{3})$ and has the derivative

$$\begin{equation*}\dot{v}(\tau )=\frac{{\dot{u}}_{\ast }(\tau )w(\tau )}{{\nu }^{2}}.\end{equation*}$$

Therefore,

$$\begin{equation*}\left[{G}^{\prime }({u}_{\ast })(v)\right](\tau )=\langle {\dot{u}}_{\ast }(\tau ),{\frac{{\dot{u}}_{\ast }(\tau )w(\tau )}{{\nu }^{2}}\rangle }_{2}=w(\tau )\frac{\langle {\dot{u}}_{\ast }(\tau ),{\dot{u}}_{\ast }{(\tau )\rangle }_{2}}{{\nu }^{2}}=w(\tau ),\end{equation*}$$

and since $w\in {H}^{1}(-1,1,\mathbb{R})$ was arbitrary, G'(u_*) is surjective. □

We can conclude now from the definition of an adjoint operator that ${G}^{\prime }{({u}_{\ast })}^{\ast }\lambda =\lambda ({G}^{\prime }({u}_{\ast }))$ where λ is a linear functional and therefore we can write the KKT condition (22) as

$$\begin{align}\hfill 0& =\underset{0}{\overset{T}{\int }}\langle y(u,t)-{y}_{m}(t),{D}_{u}y(u,t)\delta u{(t)\rangle }_{{n}_{E},2}\enspace \mathrm{d}t+\lambda (\langle \dot{u},\delta {\dot{u}\rangle }_{2})\hfill \\ \hfill & \quad +\underset{-1}{\overset{1}{\int }}{\alpha }_{1}\langle u(\tau )-{u}_{\text{ref}}(\tau ),\delta u{(\tau )\rangle }_{2}+{\alpha }_{2}\langle \ddot{u}(\tau ),\delta \ddot{u}{(\tau )\rangle }_{2}\enspace \mathrm{d}\tau \quad \forall \enspace \delta u\in {H}^{2}(-1,1,{\mathbb{R}}^{3})\hfill \\ \hfill 0& ={\Vert}\dot{u}(\tau ){{\Vert}}_{2}^{2}-{\nu }^{2}\quad \;\text{for}\;\;\text{a.e.}\;\enspace \tau \in [-1,1]\hfill \\ \hfill \lambda & \in {H}^{1}{(-1,1,\mathbb{R})}^{\ast },\hfill \end{align} \tag{ 26 }$$

with

$$\begin{align*}\hfill {D}_{u}y(u,t)(v)& =\underset{-1}{\overset{1}{\int }}\langle \nabla \omega (u(\tau )),v{(\tau )\rangle }_{2}\nu {\hat{\rho }}_{l}(\tau ,t)\enspace \mathrm{d}\tau +2\langle \nabla \omega (u(0)),v{(0)\rangle }_{2}{I}_{m}(\nu \cdot ({t}_{0}-t))\hfill \\ \hfill & \quad -\left(\langle \nabla \omega (u(-1)),v{(-1)\rangle }_{2}-\langle \nabla \omega (u(1)),v{(1)\rangle }_{2}\right){I}_{m}(\nu \cdot (1+{t}_{0}-t)).\hfill \end{align*}$$