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

In medical treatment it can be necessary to know the position of a motor unit in a muscle. Recent advances in high-density surface electromyography (EMG) measurement have opened the possibility of extracting information about single motor units. We present a mathematical approach to identify these motor units. On the base of an electrostatic forward model, we introduce an adjoint approach to efficiently simulate a surface EMG measurement and an optimal control approach to identify these motor units. We show basic results on existence of solutions and first-order optimality conditions.


Introduction
In the human body muscles are responsible for movement. These muscles consist of many muscle fibers, which are organized in so-called motor units. A motor unit is thereby the smallest controllable unit of the muscle. When using a muscle, one or more of these motor units are activated by the peripheral nervous system. This activation causes electrical signals, so-called action potentials, to propagate along the muscle fibers. These propagating action potentials create a spatially and temporally changing potential field. This potential can be measured by electrodes placed on the skin above the muscle. * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
A fundamental question in medical research and diagnosis is: what is the bioelectric source that is responsible for a specific measured potential on the skin? To answer this question, we need to find a way to identify these sources from the given measurements. Such an identification of bioelectric activity from surface measurement is needed in many fields of medicine, e.g., in measuring brain activity (EEG) or cardiac activity (ECG). Correspondingly, a lot of work has been performed to develop tools for computational assistance, see, e.g., [16] and references therein for EEG. In general terms, refinements of classical Tychonov regularization techniques are applied to the spatial problem.
The corresponding technique for measuring action potentials in muscles is called electromyography (EMG). It can be used for purpose of research (which motor-unit is responsible to which movement?) or in pre-operative planning. (Where is the location of important nerves, which should not be harmed in operations?) Similar techniques as described above have also been applied to EMG measurements [21,32,33]. These techniques yield a smooth, distributed reconstruction of sources, which is appropriate in those applications where the sources are smoothly distributed within the tissue. In these approaches, mainly spatial problems are solved, not taking into account the spatio-temporal structure of the problem.
However, bioelectrical sources in motor units are known to have special structure. They consist of a characteristic action potential that is concentrated along a line in space and moves along a fiber. Approaches that directly attempt to process the spatio-temporal information and exploit the high temporal resolution of surface EMG are less common. A notable exception is [24]. The authors consider a regularized least-squares approach for fitting the EMG signal by a linear combination of a moderate number of analytically predefined and prelocated waveforms. This enables fast computations in real time with modest accuracy.
The aim of our work is to establish a mathematically sound approach to the EMG problem that, similarly to [24], takes into account the special structure of moving action potentials in muscles. Similar to [24], we use a least-squares tracking type functional for the identification. However, we introduce refined modeling approaches to simulate a surface EMG measurement from a given source. We also represent the source to be identified more flexibly via a curve that can be chosen freely inside the muscle tissue.
Using a quasi-static approach, cf [27] the potential Φ, which is generated from the moving action potential, can be modeled for each time instant t as the solution of a Poisson equation of the form (1) In this setting the time dependent source density ρ(t), which is given through a moving action potential, is spatially concentrated on the motor unit and thus modeled as a Radon measure, concentrated on a line. Solutions to this problem in the sense of Stampacchia can be found in W 1,p (Ω) where p > d/(d − 1), Ω ⊂ R d (cf e.g., [29] and a discussion concerning uniqueness can be found in [25]). However, the direct numerical solution of (1) for all t in the time interval of interest would incur high computational cost, i.e., for each time instant an elliptic PDE on a three-dimensional domain would have to be solved. That would render numerical approaches to the identification problem too costly. We overcome this difficulty by an adjoint approach and reduce the simulation of a single measurement via (1) to the evaluation of a line integral.
Based on this simulation model, we establish a least-squares type tracking problem to identify the motor units from a surface EMG measurement. Subject of our identification problem is the trajectory of the moving action potential, represented by a parameterized curve. Unlike the least-squares approach in [24], our problem is formulated in an infinite-dimensional function space setting, which makes the analysis of the problem much more involved. Those kinds of infinite-dimensional optimization problems emerge in many application-related problems. Therefore, the prototype of this problem is studied very well, cf e.g., [18,31,36] and we will employ techniques of analysis established in the field of research to show existence of solutions and first order optimality conditions. This lays the groundwork for an optimization based numerical approach to our identification problem. A detailed elaboration of such an approach, however, will not be part of this paper but is subject to current research.

Modeling of a surface EMG measurement
Our first aim is to simulate a surface EMG measurement. Thus, for a given electric charge ρ(t) we would like to simulate the measurement of an electric surface potential y(t) ∈ R at an electrode, located on the skin.
For this simulation, we first have to derive mathematical models for the quantities of interest. First, we describe how those parts of the human body, where the motor unit lies, can be modeled mathematically. Then there is a physical model that describes the connection between the potential we want to simulate and the source that is responsible for it. We also need a model for the action potential, which acts as the source term.
As we will see, there are two problems when simulating a surface EMG measurement. First, a straightforward simulation of the measurements via solving (1) for many time instances is computationally expensive. To overcome this problem, we will introduce an adjoint approach. Second, there appear so-called end-effects, modeling artifacts, which can disturb the identification. Those end-effects are the result of violating the conservation of charge. Therefore, we add correction terms that reduce the violation.

Geometry and tissue properties
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 Ω ⊂ 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 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 μ 0 . 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 0 = 8.8e − 12 A s V −1 m −1 . The conductivity is positive and constant in these tissue type and thus the two estimates 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.
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 ∈ H 2 (−1, 1, 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 quasi-static model of the electric potential
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: 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 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: Combining these aspects, formal integration by parts yields the following symmetric bilinear form Due to (2), (3), and the presence of Robin boundary conditions with μ > 0 on ∂Ω S , this bilinear form is H 1 (Ω)-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 Since H 1 (Ω) is reflexive, we may identify H 1 (Ω) * * ∼ = H 1 (Ω) and also consider the adjoint operator A * : As already mentioned, we will model the electric charge at a time-instant t by a measure ρ(t) ∈ M( Ω). Here M( Ω) is the Banach space of Radon measures on Ω which is isomorphic to the dual C( Ω) * of the space of continuous functions by the well known Riesz representation theorem. Thus, we may introduce the following weak form: Since H 1 (Ω) → C( Ω) in our 3D-setting, we cannot write (6) as an operator equation , 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 (Ω) for some p < 3/2 can be established.
In this approach the bilinear form (5) is redefined on different spaces as . This gives rise to the following restricted pre-dual problem for some l ∈ W 1,p (Ω) * → H 1 (Ω) * : By Lax-Milgram, this problem has a solution ψ ∈ H 1 (Ω), 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 (Ω) * , which is known as 'maximal regularity', then the pre-dual operator * A p : is an isomorphism by the open mapping theorem. [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].

Remark 1. More generally, it can be shown that ψ is an element of
For simplicity we thus impose the following assumption: Assumption 2.1. The domain Ω and its subdomains Ω j are sufficiently regular, such that the operator * A p : Under this assumption and by reflexivity of Sobolev spaces, we can conclude that the adjoint of A p := ( * A p ) * 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 (Ω) → C( Ω) we can use the corresponding adjoint embedding C( Ω) * → W 1,p (Ω) * to regard the charge ρ(t) as an element of W 1,p (Ω) * for each t, and we obtain unique solvability of the operator equation: Hence, a unique electric potential Φ(t) ∈ W 1,p (Ω) 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 = * A p .

Simulated measurements by an adjoint approach
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) ∈ R from given ρ(t).
In our setting, the potential is measured with small circular electrodes on the skin as follows 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 (Ω) * . Let ρ(t) ∈ W 1,p (Ω) * and denote by Φ(t) ∈ W 1,p (Ω) the solution of Now consider the solution ω ∈ W 1,p (Ω) of the adjoint problem: which corresponds to the following Poisson problem in weak form: Then we compute easily That means we can compute the potential at an electrode efficiently by evaluating where ω ∈ W 1,p (Ω) → C( Ω) 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 ω 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 ω| Ω 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).
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 ω| Ω M to Ω M . We observe that ω| Ω M satisfies a homogenous Laplace equation in weak form on Ω M subject to Dirichlet boundary conditions on ∂Ω M , given simply by the condition that ω| ∂Ω 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].

Moving action potentials
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.
Biomedical modeling of action potentials starts with the following function in terms of a reference parameter z ∈ 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 which corresponds to the principle of conservation of charge in the body. Up to now, the action potential is defined as a function on 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 |u(τ )| ≡ ν. That means we can identify each point on the curve u with some z ∈ R via the arc length 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 ∈ R + and points on the 'left' side with some z ∈ 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 0 − t) and set 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: 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 R vanishes. However, our definition of y(t) involves only an integral over a bounded subset of R and the corresponding total charge is given by (taking into account the substitution of variables formula): 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 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: Observe that these charges are 0, if t 1 + t 0 , 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: 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).

Numerical simulation of a surface EMG measurement
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.
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 T of Ω consisting of 414 195 tetrahedra. On T we used continuous piecewise quadratic ansatz functions to discretize W 1,p (Ω) and W 1,p (Ω) by A Galerkin method, applied to the adjoint problem (8) leads to the discrete problem 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 ∈ 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 pointx on the line segment [x i , x i+1 ] for which a neighborhood search is successful:x ∈K. Then we repeat the whole procedure with (x i , K i ) replaced by (x,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.  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.

Identification of a motor unit from measurements
In this section we will specify the identification problem which we want to solve. It is the inverse problem to the forward problem described in section 2 in the following sense: up to now, the motor unit was modeled as the trace of a given curve u and we derived a model for the simulation of the measurements y i (u, t) at n E electrodes via (16). From now on, we assume that measurements y m,i (t) are available and we are looking for a curve u, such that the corresponding simulated values y i (u, t) and the measurements y m,i (t) fit well. Collecting all these measurements and simulations in the vectors y(u, t), y m (t) ∈ R n E and using the standard Euclidean norm · n E ,2 on R n E , this leads to the following least-squares type tracking term n E ,2 dt : the components of y(u, t) are defined via (16).
A rough guess of the location of the motor unit u can be made by inspecting the given measurement y m and the subdomains of Ω. Thus, we can choose a reference trajectory u ref (e.g. a piecewise linear curve that connects the estimated location of the neuromuscular junction and the end-plates) a prior and add the following regularization term to the problem (where · 2 is the standard Euclidean norm on R 3 ): Additionally, since motor units are smooth in healthy tissue we add a second regularization term, which is given through This term also yields the necessary compactness the show existence of optimal solutions. Finally, we will add a constraint which ensures that the signal passes the motor unit with constant speed ν > 0, as assumed in section 2.4. Therefore, we define the constraint function and demand that [G(u)] (τ ) = 0 for almost every τ ∈ [−1, 1]. We also demand that the solution is located in the muscle tissue Ω M . Combining those two constraints, we get the following admissible set We note that this definition of the admissible set makes sense, due to the fact that H 2 (−1, 1, R 3 ) is embedded in C 1 (−1, 1, R 3 ). Collecting everything we get the optimization problem Alternatively, we can write this problem as an unconstrained problem by adding an indicator function, such that we get Now that we have derived an identification problem, we want to know if the problem has a solution.
Remark 2. Some straightforward extensions of this identification problem are conceivable: for example we may include the speed ν in the set of variables, to be identified, and similarly the scaling parameter a. For simplicity of presentation we assume these parameters to be given.

Existence of a solution
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.
Thus U 1 is well defined and not empty.
Let now {u k } ⊂ U 1 be a weak convergent sequence with limit u. Since the embedding E 1 is compact there exist a subsequence u k l such that E 1 u k l → E 1 u in C (−1, 1, R). Furthermore, there exist an another subsequence E 1 u k l i that converges pointwise to E 1 u for all τ ∈ [−1, 1].
Since Ω is closed and we can conclude from E 1 u k l i (τ ) ∈ Ω that E 1 u(τ ) ∈ Ω and thus also u(τ ). This shows that U 1 is weakly closed.
As above there exits a compact embedding E 2 : H 1 (−1, 1, R 3 ) → W 1,4 (−1, 1, R 3 ) with E 2u =u. We can conclude from Hölder's inequality that u,u 2 ∈ H 1 (−1, 1, R) and thus is well defined. It is well known that, as a continuous bilinear form, G is Fréchet differentiable with derivative Let now {u k } ⊂ U 2 be a weakly convergent sequence with limit u. As above there exist a subsequence u k l such that Eu k l → Eu in L 4 (−1, 1, R 3 ). Since G is differentiable it is also continuous. This means that

Remark 3.
The H 2 -regularization term is essential for weak closedness of U 2 . It is not hard to construct a zig-zagging sequence of trajectories u n with |u n | = ν a.e. (and thus bounded in W 1,∞ ), such that the weak limitū violates |u| = ν. From a computational point of view, such a regularization term does not impose severe difficulties, since u can be easily discretized in an H 2 conformal fashion by a piecewise polynomial spline which is globally C 1 . Proof. First we note that, as a sum, J is continuous, Fréchet differentiable and weakly lower semi-continuous, if J 1 , J 2 and J 3 are continuous, Fréchet differentiable and weakly lower semi-continuous.
To show the three properties for J 1 , we define for fixed t and k the mapping From lemma 2.2 we now that ω k : Ω M → 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 : [5, theorem 8.8]. Therefore, the superposition operator is well defined and Fréchet differentiable, see [2, theorems 6.3 and 6.7]. The Fréchet derivative is given by By the same argumentation the correction terms ω k (u(·))I m (·) are well defined and Fréchet differentiable. Thus, is well defined and Fréchet differentiable with derivative The chain rule indicates that J 1 (u) = 1 is differentiable, it is also continuous and thus lower semi-continuous in C(−1, 1, R 3 ). By the compact embedding E : H 2 (−1, 1, R 3 ) → C(−1, 1, R 3 ) we conclude that J 1 is weakly lower semi-continuous in H 2 (−1, 1, R 3 ). J 2 and J 3 are obviously convex quadratic bilinear forms and it is well known that they are Fréchet differentiable from L 2 into itself and weakly lower semi-continuous, cf e.g. [9]. The derivatives of J 2 and J 3 are then given through 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

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 → R is weakly lower semi-continuous if the level sets N α f = {x ∈ X| f (x) α} are weakly closed for all α ∈ R, see [7, theorem 7.4.11]. For the indicator function the level sets are given through 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 2 , J 3 or ι U 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: 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, 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 inf F F(u * ) lim inf k→∞ F(u n k ) = inf F, which shows that the limit point u * is a minimizer of F.

First-order optimality conditions
Next we derive first-order optimality conditions for the problem. This is complicated by the geometric constraint u ⊂ Ω 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 2 of the objective functional is penalizes by J 2 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: As usual we eliminate the equality constraint with the help of a Lagrange multiplier, which lead to the following result: Proof. First, we recall that λ ∈ H 1 (−1, 1, R) * is called Lagrange multiplier, if λ (G(u * )) = 0, 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 = {0} ⊂ H 1 (−1, 1, R) and from the definition of the dual cone we can conclude that K + = H 1 (−1, 1, R) * , see [31,36]. [G(u * )] (τ ) = 0 implies that condition (24) is always fulfilled and thus we can replace it by [G(u * )] (τ ) = 0. Furthermore, C(u * ) is the conical hull of H 2 (−1, 1, R 3 ), which is H 2 (−1, 1, R 3 ) and therefore the dual cone C(u * ) + = {0}. Thus, condition (25) simplifies to J (u * ) − G (u * ) * λ = 0.
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 G (u * )C(u * ) + K(−G(u * )) = H 1 (−1, 1, R).

Conclusion and outlook
We were able to convert an identification problem in EMG into a tractable optimization problem in function space. A decisive step was the use of an adjoint approach to enable an efficient simulation of the measured signal. The derived first-order optimality conditions (26) serve as basis for a computational approach, which we plan to elaborate in a forthcoming paper. In this approach the weighing functions ω i are precomputed by finite elements and the curve u is represented by a C 1 piecewise polynomial spline. An optimization algorithm that makes use of second order models of the objective can be applied to the resulting problem.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.