Model Order Reduction with Dynamically Transformed Modes for Electrophysiological Simulations

: The numerical simulation of neuromuscular processes in skeletal muscles typically imposes high computational costs, resulting in the need for techniques to reduce these costs. One submodel in these types of multi-physics multiscale simulations is the electrophysiological model, describing the propagation of an action potential (AP) in a muscle fiber. Involved in the simulation of the propagation are two traveling waves, resulting in difficulties for classical model order reduction (MOR) techniques based on linear subspaces. Instead, we apply the nonlinear MOR method shifted Proper Orthogonal Decomposition (sPOD) to one state variable only in order to construct dynamically transformed reduced basis functions depending on time-dependent paths spanning the adaptive reduced ansatz space. Our numerical experiments demonstrate that the constructed reduced ansatz space can accurately capture the dynamics of two fully separated wavefronts and reduce the degrees of freedom of the whole simulation. However, it cannot represent overlapping wave parts of the hyperpolarization. The constructed reduced order model outperforms the high-dimensional full order model in terms of the computational costs while the accuracy is maintained and reaches speedup factors between 2 and 73 depending on the time discretization.


Introduction
A majority of tasks in human life can only be accomplished because of the complex muscular processes and interactions of and between the approximately 220 individual muscles [17].At the core of all movements are voluntary contractions of skeletal muscles, allowing for a large variety of different movements, covering fast, persistent, finely controlled, but also exhausting actions.This diversity is based on complex interactions between the musculoskeletal and the neurological system.
Gaining insights into these systems is of immense importance, for example, in medicine and healthcare.Most incurable genetic disorders like myotonia con-genita stem from a malfunction in the electrical processes regulating muscle contractions [23,31].A better understanding of the mechanisms of the neuromuscular system can improve treatments and classification, also for other muscular diseases and in post-stroke rehabilitation using patient-specific motorized prostheses [24,33].In the field of engineering, the development of bioinspired robots and artificial muscles, mimicking the functionality of the human skeletal muscle, requires deep knowledge about the involved neurological processes [21,30].Increasing the functionality and usability of such biomimetic systems can have advantages in industrial robotics and the development of exoskeletons [2,3].
Traditionally and empirically oriented approaches in clinical studies have reached their limits regarding the level of detail and the time requirements [10].Instead, individualized, biophysics-based simulations can be used to gain detailed insights into non-measurable processes.The numerical treatment of the underlying coupled, nonlinear partial differential equations (PDEs) necessitates large-scale computation facilities to manage the extensive computational costs required for high fidelity solutions.These computational costs are not exclusively a problem in research because of long computation times and high energy consumption.In envisioned real-life applications, in medicine for example, patient-specific simulations require low runtimes of the detailed simulations to supply fast and exact diagnostics, for example in adaptive surgery settings.To allow for execution on resource-constrained devices and the treatment of the envisioned scenarios, the computational costs and the complexity of such simulations must be reduced.

Outline and Contributions
In this work, the simulation of neuromuscular processes is based on a multiphysics, multiscale model capturing the processes in the individual muscle fibers to the contraction of the whole muscle, see Section 2.1.Since the number of fibers comprising a skeletal muscle is typically high, the simulations of the individual fibers are the dominating factor in the overall computation.Simulating each fiber involves the electrophysiological model capturing the propagation of the cell membrane potential, defined as action potentials (APs).Reducing the computational complexity of this complex electrophysiological model is of utmost interest.
One promising approach to reduce the computational costs is model order reduction (MOR).Generally, in MOR, the full order model (FOM) is replaced by a suitable reduced order model (ROM) with lower computa-tional complexity while ensuring the accuracy of the results as much as possible [4].
Recent works in the field of MOR for electrophysiological simulations involve Proper Orthogonal Decomposition (POD) and Deep Learning (DL) and show the potential of these techniques for this kind of application.In contrast to this work, these approaches were applied to cardiac electrophysiology, the simulation of the electrical activity of the heart, and not to single skeletal muscle fibers [9,13,14,16].
This paper focuses on the derivation of a ROM with transformed modes, based on shifted Proper Orthogonal Decomposition (sPOD) for an electrophysiological simulation of skeletal muscle fibers.In detail, the contributions of this work are: • The application of sPOD to construct a reduced ansatz space capturing the wave propagations.• The use of this reduced ansatz space in a finite element approach to obtain the ROM.• The evaluation of the reduced ansatz space and the performance of the reduced simulation.
In doing so, this work shows that the application of sPOD in the context of MOR for electrophysiological simulations is a valid approach.Even by applying the MOR approach to one state variable only, speedup factors between 2 and 73 can be reached with maintained accuracy.
The remainder of this paper is structured as follows.Section 2 deals with the necessary background of the simulation model, the construction of the FOM via the finite element method, and its simulation.Section 3 introduces the MOR approach, followed by the details of the reduced ansatz space and the resulting ROM in Section 3.4.Performance and accuracy of this MOR approach are evaluated in Section 4. Finally, Section 5 summarizes the results and gives an outlook on future work.

Full Order Model
Before we present the technical details of our MOR approach, we first give some background to the original full order simulation of the propagation of the APs.
The APs arrive at the neuromuscular junction located in the center of each fiber and propagate along the fiber to both ends.This propagation is triggered by biochemical processes on the subcellular scale involving ion movement through voltage-dependent ion channels.To capture these complex processes we use the overall model described in [10,27] and extended in [26].
The model is based on two submodels, the 1D fiber activation model and the 0D subcellular model.The first one captures the propagation of the AP along a 1D fiber, whereas the second one describes the chemical reaction to this potential in the fiber cell.Hence, the 1D fiber activation model and the 0D subcellular model are strongly coupled, and we combine them in the electrophysiological model.The 3D mechanics of the actual muscle movement are not covered in this paper, and we refer to [26].

Electrophysiological Model
The mathematical definition of the electrophysiological model is based on [29] and described in [26].Its derivation is mainly based on the electrical coupling between intra-and extracellular spaces of the muscle fiber and the Hodgkin-Huxley cell model, proposed in [19].
For each bounded fiber domain Ω = [0, ℓ] with the muscle fiber length ℓ and simulation time 0 ≤ T 0 < T , the 1D model captures the transmembrane potential u : Ω × T 0 , T → R along the fiber and is a second-order PDE of the form The first term in the 1D model describes the diffusion of the electrical potential along the fiber and the second one accounts for the transport of the potential, driven by the ion movements.This potential-dependent current of the ion movements I ion can be described by the model of Hodgkin-Huxley [19].The involved constants are the fiber's surface-to-volume ratio A m , the capacitance of the muscle fiber membrane C m , and the effective conductivity σ eff .
The Hodgkin-Huxley model describes the movement of the k > 0 ions in terms of ion-specific conductivities g i (y) and constant Nernst potentials E i : These conductivities depend nonlinearly on the activation probabilities of the ion channels y : Ω× T 0 , T → R k and describe the pass-through of the ions through the membrane.Typically, k ≈ 50 activation probabilities are involved in capturing the complex processes in our muscle fibers.These probabilities are defined locally on a subcellular scale, more precisely at approximately 100+ positions on the muscle fiber cell.Depending on the transmembrane voltage, the activation probabilities change over time.This process is captured by the 0D subcellular model at every point on the fiber, formulated in the Hodgkin-Huxley model as well Prescribing initial and boundary conditions as follows ∇u = 0, on ∂Ω, yields the final initial boundary value problem (IBVP) of the electrophysiological model.Formally, we require

Finite Element Approach
The original full order simulation on the 1D fiber scale is based on a finite element method and is used as a reference solution of the electrophysiological model problem from Section 2.1.Constructing this finite element method results in the semi-discrete FOM, yielding a finite-dimensional approximation to the weak solution of the model problem.
For a fixed time t ∈ T 0 , T , we seek the spatial weak solution u(t ) in the Sobolev space Y = H 1 (Ω), more precisely with (u(t ), y(t )) ∈ Y × Y k the weak forms of (1) and ( 2) must hold for all v ∈ Y and w ∈ Y k .To simplify the notation, we take Y as the ansatz space for y(t ) intentionally, even though L 2 (Ω) would suffice.We cover the spatial domain Ω by a triangulation T h as a subdivision of equidistant subintervals with length h s and use nodal basis functions ψ 1 , ..., ψ n to construct the finite dimensional ansatz space V h where we restrict ourselves to the Lagrange P 1 (T h ) space for simplicity in the following.Performing a Galerkin approximation in space for u h ∈ V h with u h = n i =1 ūi ψ i and analogously y h ∈ V k h , we arrive at the semi-discrete formulation with matrices and vectors The identity matrix is denoted as I k ∈ R k×k and ⊗ corresponds to the Kronecker product.The nonlinear terms from (1) and ( 2) are represented in the ansatz space as the maps f : To solve for the time-dependent solution ū and ȳ, we apply an explicit Euler time-stepping scheme for simplicity.
For this, we partition the time interval [T 0 , T ] into m equidistant subintervals with length h t .

Model Order Reduction
The simulation of the FOM is computationally expensive due to the coupling between the two complex submodels and a large number of discretization nodes both in space and time.In this section, we present a MOR approach with the general aim to reduce the computational costs and allow for rapidly computable approximations, while still achieving a satisfactory accuracy of the FOM, building on ideas from [4].AP wave propagation for a muscle fiber of length 12 cm, originating in the center of the fiber.FOM simulation with the parameters from Table 1 and h t = 0.001.
Our approach is based on the construction of a problemspecific, time-dependent, low-dimensional subspace of the full dimensional ansatz space and using this ansatz space to derive the ROM.More precisely the reduced ansatz space should be able to capture the dynamics of the solution, in our case the two traveling wavefronts in Figure 1.We focus on constructing and applying our MOR approach for the transmembrane voltage u for simplicity, not for the ion channel probabilities y.

Overview of MOR Approach
Since our problem is dominated by transport, we apply the sPOD introduced in [35,36] and formalized in [5,6,8].To capture the transport terms, sPOD allows the construction of multiple reduced ansatz spaces traveling in the spatial domain.But instead of moving the discretization nodes like in Moving Finite Element Methods [15,28], the reduced basis functions {ϕ i } r i =1 evolve over time themselves and the discretization is fixed.This can be achieved by shifting the modes over time, introducing an additional variable p j (t ) for the shift amount.From the full order simulation in Figure 1 we observe that the wavefront of the transmembrane voltage travels at an approximately constant speed through the spatial domain.Following this observation, we only account for constant shift velocities of the basis functions and construct a constant-speed shifted reduced ansatz space, in contrast to [5,8] where also non-constant shift velocities are allowed.
We use sPOD [8,34], an adapted version of POD [18,20,22,32], to obtain the reduced basis functions from the results of the full order simulation, explicitly accounting for the shift variables p j .To obtain the ROM, we use these reduced basis functions instead while following the same Galerkin approach in comparison to the FOM.

Adaptive Reduced Ansatz Space
With the general goal to construct a ROM, we aim to approximate the solution u(t ) ∈ Y by the reduced solution û(t ) ∈ R in some reduced space R ⊂ Y .We construct the basis functions {ϕ j } r j =1 of the reduced space R as elements of the finite dimensional ansatz space V h of the FOM to obtain continuous reduced basis functions.Since they should capture the transport quantities, we evaluate them in a shifted and time-dependent position based on unknown paths p j : T 0 , T → R [8].
Since the path functions describe the evolution of the ansatz space and thus should capture the constant velocities we can assume that the paths p ∈ H 1 (T 0 , T ; R r ) satisfy With these preparations, we can define our modified reduced ansatz space.Let T (p j (t )) be a suitable transformation operator, compare [8], with linear shifts p j (t ) = p j t with the path coefficients p ∈ R r for all j = 1, ... , r .We define the time-shifted reduced basis functions {φ j } r j =1 through and use the periodic extension proposed in [5,Ex. 4.3 and 5.2] to keep the shifted basis functions ϕ 1 , ..., ϕ r well-defined.With these time-shifted basis functions we span the constant-speed shifted ansatz space Vh (t ) with Vh (t ) = span φ 1 (t ), ... , φ r (t ) , and r = dim Vh (t ) < dimV h = n.This ansatz space travels in the spatial space Ω depending on time.

Constant-Speed Shifted Proper Orthogonal Decomposition
To use the constant-speed shifted reduced ansatz space explicitly in our ROM, we need to determine the basis functions φ j .Therefore, we follow the general idea of snapshot decomposition to solve a similar minimization problem as presented in [8].The differences are hidden in the details: Since we assume linear paths, it suffices to optimize over the constant path velocities p instead of the paths themselves.However, the constant shift velocities p must be optimized nevertheless.
The available snapshot data in our case are function samples from the FOM simulation and denoted as u ∈ L 2 (T 0 , T ; Y ).In addition, let L := L 2 (T 0 , T ; R r )×Y r ×R r be the space of the optimization variables.Since the overall goal is to represent the solution u in the reduced space, the error of the representation should be as small as possible.To this end, we define the cost functional J : L → R as the mean squared distance to the function sample u .
The optimization problem then reads: holds.We solve this optimization problem using a quasi-Newton method for limited memory, the L-BFGS-B algorithm [12,25].Similar to [8], we provide the explicit discretized cost function and the corresponding gradient by using the same discretization from the FOM.

Reduced Finite Element Method
Solving the discrete optimization problem from Section 3.3, we obtain the coefficients pi and φ to construct the shifted ansatz space Vh .Similar to the derivation of the FOM, we insert the reduced basis representation of the Galerkin approximations ûh = r i =1 ûi φ i and y h = n i =1 ȳi ψ i into the weak formulation of the original problem.Here V h = P 1 (T h ) is the full-dimensional discrete ansatz space of Y with nodal basis ψ 1 , ..., ψ n .The resulting semi-discrete ROM is thus given by with the time-dependent matrices and vectors To assemble the reduced matrices, we need to compute the L 2 (Ω)-inner products with the time-shifted basis functions φ(t ).The major difference in comparison to the FOM from Section 2.2 is the time dependency of the basis functions.Since the paths account for the constant velocity shift and are thus given by the reduced space, they do not occur as unknowns of the ROM thus the paths are not recalculated in the online phase.

Computation of the ROM
Solving the semi-discrete ROM from Section 3.4 involves the application of the explicit Euler scheme, similar to the FOM.This time-stepping scheme depends on the model matrices and for the FOM it is sufficient to assemble these model matrices once due to their time independence.However, the model matrices for the ROM depend on the shifted reduced basis functions and thus on the time-dependent paths.This requires their reassembly in every time step, which can be computationally expensive.
In our implementation, we evaluate the terms Fh and Ĝh using the nonlinear terms F h and G h , defined on the full-dimensional ansatz space V h .To do so, we have to project the transmembrane voltage from the constant-speed shifted ansatz space ûh ∈ Vh to the full order ansatz space u h ∈ V h .We call this up-projection and use an L 2 projection to define it.Starting at the condition 〈u h − ûh , ψ i 〉 L 2 = 0 for all i = 1, ..., n we can obtain the coefficients ū of the basis representation in the full order space V h based on the coefficients û of the basis representation in the reduced space Vh .This transformation is defined as with the projection matrix Ph Based on these coefficients we can easily represent a function ûh = n i =1 ûi φ i ∈ Vh in the space V h through the basis representation u h = n i =1 ūi ψ i .Following the same procedure, we can also construct the so-called down-projection mapping from V h to Vh with Ph = P T h to obtain the evaluations Fh and Ĝh .These additional steps must be taken into account in the algorithm when solving the ROM.
Typically, the evaluation of the ROM still scales with the dimension of the FOM and for that reason hyperreduction techniques like discrete empirical interpolation method (DEIM) are commonly used to reduce the computational complexity.However, in our setting, we solely focus on the sPOD without applying hyperreduction techniques and note that the performance of the ROM can be further improved by using these methods.

Numerical Experiments
In this section, we present the application of the MOR scheme defined in Section 3 to the electrophysiological model.All computations are performed on an 11th Gen Intel® Core™ i7-11700 @ 2.50GHz with 8 cores and a total DDR4 RAM capacity of 16 GB @ 3200 MHz.The implementation in Python is strictly sequential and the sparsity of the matrices is exploited by using CSR matrices.
The accuracy of the presented methods, especially of the constructed reduced basis in Section 4.2 and the resulting ROM in Section 4.3, is evaluated based on the relative error in the L 2 -norm and the maximal error in the maximum norm.We define the global errors between a function z ∈ L 2 (T 0 , T ; L 2 (Ω)) with ∥z∥ L 2 (T 0 ,T ;L 2 (Ω)) ̸ = 0 and its approximation as follows We refer to the local absolute error as e abs and define it as e abs (t , for almost all (t , x) ∈ [T 0 , T ] × Ω.In addition, we denote the arithmetic mean of the absolute error over the entire spatial and time domain as ēabs , i.e.
If not stated otherwise, all the following errors refer to the transmembrane voltage state z = u only, where u is the solution from Section 4.1.To compare this function with approximations in the reduced setting, we project the reduced functions to the full dimensional space where û are the amplitudes in the basis representation of the reduced space and ψ is the vector containing the nodal basis functions of V h .Furthermore, we differentiate the errors in offline errors and online errors.The offline errors are the errors of the reduced ansatz space based on the optimal amplitudes resulting from the optimization problem (4), which correspond to the orthogonal projection onto the (time-dependent) reduced ansatz space.These offline errors are used to evaluate the quality of the results from the sPOD.We denote the online errors as the errors based on the solution of the ROM as amplitudes.This separation allows analyzing if the errors stem from the ROM or the reduced ansatz space.

FEM Simulation
The model parameters used in all our simulations are taken from [26] and listed in Table 1.In all the simulations we use k = 3 as the number of voltage states and for the simulation of the FOM we use a time-step ith µ = ℓ/2, ς 2 = 0.02 cm 2 and V max = 35 mV for all x ∈ Ω.Furthermore, we set the ion channel states to their equilibrium probabilities given in Table 1 with y 0 = P Na + , P K + , P L T for all x ∈ Ω.
The resulting voltage snapshots are visualized in Figure 1 and show the AP defined by the voltage as two wavefronts traveling through the spatial domain in both directions from its origin.The solution for all states can be seen in Figure 2, where the propagation in time is visualized by the color gradient ranging from red (t = 0) to green (t = T ).

Quality of Reduced Ansatz Space
The reduced ansatz space should capture the propagation of the voltage state in the spatial domain, and thus we only use the snapshots for the voltage to construct the reduced basis functions.To determine these reduced basis functions we numerically minimize the cost functional (3).This subsection corresponds to the offline phase and hence all errors are offline errors.The optimization problem ( 4) is carried out by SciPy's minimize function using the L-BFGS-B algorithm [37].
Based on the numerical experiments in [7], we expect it to be challenging to achieve a good approximation of the initialization of the wave, also called activation, when using a sPOD-based ROM.For this reason, we focus on determining the reduced basis functions for the cases where the waves are separated and extract the initial values and the snapshot data from the full order simulation results for different start times T 0 .
Typically, an AP consists of a narrow peak at the wavefront and a wide hyperpolarization part afterwards, characteristic for the decrease of the potential below the resting potential [1].We furthermore differentiate between two different cases of the snapshot data.In the first one with T 0 = 5 ms only the wavefronts are separated but the hyperpolarization parts overlap.In the second one with T 0 = 18 ms the waves are almost completely separated, and the hyperpolarization parts only slightly overlap.
To reduce the offline computation time for sPOD, we do not use the whole time domain of the snapshot data for the optimization problem.Instead, we use the snapshot data until the fourth millisecond, respectively for the two cases T 0 , T 0 + 4 .
Generally, we expect the cost functional of the optimization problem of sPOD to be highly non-convex.Hence, a good choice of initial values is crucial for the success of sPOD.Since we assume a constant propagation speed, we initialize the paths by linear regression over the whole time domain separately for both waves, similar to [7, Algorithm 1].

Completely Separated Waves
We start with the less challenging scenario based on the snapshot data for T 0 = 18 ms with completely separated waves and two basis functions.For the initialization, we split the snapshot data at t = T 0 in the center of the domain and use each half separately for different basis functions with the other half evaluated to zero.We take these values as the amplitudes in the linear combination of the basis functions ψ 1 , ..., ψ n to construct the initial basis functions ϕ 1 , ..., ϕ r .In accordance with the theoretical construction in Section 3, these basis functions are continuous over the whole domain Ω.
We obtain a relative error of e rel = 0.405 and a maximal absolute error of e ∞ = 102.646.From comparing the coefficients of the initial paths p (0)  i and the optimized paths p i , p(0) 1 = −0.245,p(0) 2 = 0.245, p1 = −0.187,p2 = 0.187, we observe that the wave is inhibited from traveling through the spatial domain to remedy the hole in the center resulting from shifting the basis functions in opposite directions to the sides.This is a general problem because the resting potential V eql is not zero, consequently, the chosen initialization is  not optimal for this case but shows an underlying problem.Shifting basis functions away from one another results in an area in the center where the basis functions evaluate to zero by construction, compare Figure 3(a).
To deal with this problem, we define the so-called equilibrium transformation that transforms the absolute voltage u(t , x) to the relative voltage u rel (t , x), such that u rel (t , x) evaluates to zero in the equilibrium state.This transformation can easily be carried out by subtracting the resting potential from the data This transformation also has to be applied during the projection between the reduced and original space, and we continue under the assumption that the equilibrium transformation is applied in any necessary case.
With the transformed setting we observe a significantly smaller relative error of e rel = 0.0027 and a maximal absolute error of e ∞ = 1.2607, see additionally Table 2.A further evaluation using the absolute errors is visualized in Figure 4.It shows that the wavefronts are captured accurately; however, a small amount of the overlapping wave parts in the center cannot be represented.
We increase the number of basis functions to r = 4 and use the two already optimized basis functions as initialization and initialize the two additional basis functions such that they evaluate to zero with a non-zero gradient.However, this leads to no improvement due to the non-convergent behavior of the optimizer and results in unchanged basis functions.In order to increase the accuracy by adding additional functions, investi-gating different approaches to the optimization problem might be necessary.
Since the cost function of the optimization problem is expected to be highly non-convex, it is likely that the optimization algorithm starts close to a local minimum by initializing the basis functions with the snapshot.To analyze this aspect, we initialize r = 2 basis functions such that they constantly evaluate to zero in the whole domain.These experiments result in the relative and absolute errors listed in Table 2 as the undivided case.To analyze this problem further, we follow the approach of separately constructing basis functions for the divided domain mentioned in [7].Therefore, we divide the snapshot data along its center axis, set the corresponding other half to the equilibrium state, construct the reduced space with one basis function r = 1 each for the left and right wave separately and merge the reduced spaces to obtain a merged space with r = 2.In that way, we obtain much smaller relative and maximal absolute errors as listed in Table 2 as well.We observe that the relative errors are similar to the results from the snapshot data initialization for r = 2 in Table 2.However, using the snapshot data initialization for r = 2 yields smaller maximal absolute errors.

Conclusions
Transforming the snapshots to the equilibrium state allows to remedy the numerical impact of shifting basis functions in opposite directions, leading to zero evaluations of the reduced ansatz space in the cen-ter of the domain.Constructing the reduced space for both wavefronts at once with zero-initialized reduced basis functions leads to high relative and maximal absolute errors.In addition, by splitting up the wavefronts and separately constructing reduced spaces for them with zero initialization, the wavefronts can be accurately captured.We conclude that with the sPOD and the chosen optimizer, the assignment of the reduced basis functions to the two different wavefronts strongly depends on the initialization.The optimizer could converge to a local minimum due to the non-convexity of the cost functional, causing this effect.

Partly Overlapping Wavefronts
We continue with the more challenging scenario of the partly overlapping wavefronts based on the snapshot data from T 0 = 5 ms.For an initialization with the snapshot data, we obtain larger relative and maximal absolute errors, listed in Table 3, compared to the scenario of the completely separated wavefronts.Figure 6 shows that these larger errors mostly stem from the center part and not from the wavefront.In addition, we can see that the wavefront is captured accurately by comparing the reconstructed waves with the original waves, visualized in Figure 7.However, by zooming into the center, we observe that the hyperpolarization part of the AP is almost completely cut off and cannot be represented by the reduced space.This behavior is similar to the reduced space for the completely separated wavefronts with snapshot data initialization, where the center part is the cause of the errors but has less impact.Following that, we can deduce that with an increase in the overlap of the hyperpolarization parts of the two waves, the error in the center increases as well and information about the hyperpolarization is lost.
For the zero initialization, we obtain the relative and absolute errors listed in Table 3. Splitting the snapshot data and constructing two reduced spaces reduces both the relative and maximal absolute error.Constructing both reduced ansatz spaces separately leads to a smaller relative error compared to the snapshot data initialization.Looking at the shape of the reduced basis functions in Figure 8 we notice that this decrease in the errors most likely stems from a better approximation of the hyperpolarization part.The initialization with snapshot data most likely got stuck in a local minimum.

Conclusions
We conclude that in this scenario, using the split snapshot as initialization allows us to capture the non-overlapping wavefront, similar to the scenario of completely separated wavefronts.In addition, the reduced space cannot reproduce the overlapping parts of the wave, especially the hyperpolarization part can-not be represented in the reduced space.Using a snapshot initialization is the best option if we are only interested in capturing the wavefronts.Thus, based on the presented approach, we can only approximate the full waves accurately after they are completely separated.

Simulation with the Reduced Order Model
With the constructed reduced space for the two scenarios, we can now evaluate the accuracy and efficiency of the ROM.To this end, we distinguish between the two scenarios from Section 4.2 with the different starting times T 0 = 18 ms and T 0 = 5 ms and use the reduced ansatz space with r = 2 reduced basis functions constructed with snapshot initialization.In the following, we determine the accuracy of the ROM in comparison to the FOM and compare these online errors with the offline errors in Table 4. Additionally, we use the computation time of the FOM simulations to calculate speedup factors, indicating the efficiency of the ROM.We round these speedup factors to the nearest integer and use the computation times of the FOM listed in Table 4.

Completely Separated Waves
We start with the first scenario (T 0 = 18 ms) for completely separated waves.We simulate the ROM for five time-step widths h t ∈ {0.001, 0.005, 0.01, 0.05, 0.1} and immediately observe the instability of the explicit Euler time-stepping scheme for a time-step width of h t = 0.1.
From the results in Table 5, we observe a slight increase in the relative errors while the computation time is drastically reduced.Even if the simulation of the ROM uses the same time-step width of the FOM, the simulation time is two times faster.
Exploiting the possibility of increasing the time-step width up to 0.05 ms, the computation time can be reduced by a factor of 73, with online errors of the same order of magnitude as the offline errors, compare Table 4. Hence, the error from the reduced ansatz space dominates the error in the simulation.We only have to accept an increase of approximately 10 % in the maximal absolute errors and an increase between 3 % and  To further investigate these errors, we focus on the solution of the ROM for the time-step width h t = 0.001 ms.From Figure 9, we observe that the center, and thus the end of the hyperpolarization part, has a large impact on the errors even for the small time-step width and can not be represented accurately.However, since we cannot expect that the ROM can capture dynamics that the reduced space cannot represent, this effect was expected in the online phase as well since it has occurred in the offline phase.Comparing the voltage wavefronts of the ROM solution and the original wavefront in Figure 10, we note that the wavefront is captured accurately.

Partly Overlapping Waves
We perform the same experiments for the more challenging part with partly overlapping waves and obtain the results listed in Table 6.Comparing the online and offline errors, we observe that they have the same order of magnitude but with 50 % larger relative online error and an approximately three times larger maximal absolute online error.In addition, we can see in Figure 11

Ion Channel Probability States
Since the voltage state directly affects the ion channel probability states, we analyze the impact of the errors in the voltage state on these other states as well.Therefore, we exemplarily use the ROM for completely separated wavefronts simulated with h t = 0.001.From the absolute error snapshots for these states, visualized in Figure 12, we observe a small error in the center for all states y 1 , y 2 , y 3 .In addition, a small absolute error at the wavefronts for the sodium ion channel probability y 1 is present.The error in the center is caused by the lack of   expressiveness of the reduced ansatz space for the end of the hyperpolarization in the center.Depending on the time step size these errors increase but the shape of the front profile of the traveling waves is maintained.

Discussion
We have applied the sPOD, which is based on dynamically transformed basis functions, to the electrophysiological simulation of the propagation of APs in muscle fibers.The reduced approximation ansatz has been based on a linear combination of transformed basis functions with constant transport velocity.The propagation of the AP through the spatial domain originates from the center of the muscle fiber, resulting in overlapping wavefronts in the activation area.The ROM with the constructed reduced ansatz space can capture the wavefront accurately with a small relative error.However, the complete wave, especially the large hyperpolarization part, can only be captured for the completely separated waves.
We have observed a drastic reduction in the simulation time compared to the FOM computation.One important aspect is the possibility of using larger timestep widths for the time-stepping scheme resulting in fewer iterations.In that way, we have shown that applying MOR techniques in biophysical simulations for medical applications is indeed a valid approach to drastically reduce the computation time while still achieving a satisfactory accuracy of the simulation results.
Despite the achieved performance under mostly preserved accuracy, some difficulties still need to be addressed in the future.First, the activation of the AP in the center of the fiber must be integrated into the ROM.One promising approach follows [7], combining the presented ROM for the wavefront with another ap-proach for the activation in the center.Another critical topic is the hyperpolarization part of the wave.Biologically, this part is essential for activating another AP and thus should be captured by the ROM.Since the hyperpolarization parts make up a large percentage of the complete wave and are nearly separated after 18 ms, using different reduced basis functions for the wavefronts and the hyperpolarization parts could increase the approximability.To enable predictions in the online phase using parameter settings that have not been used in the offline phase, the ROM has to be extended to paths being treated as time-dependent unknowns of the ROM, instead of being fixed.In that way, the propagation depends on the other states of the ROM, similar to the wildland fire simulation in [7].In addition, extending the ROM to the ion channel probabilities in the electrophysiological model can also potentially increase the simulation performance further.Despite that the ion channel probabilities behave similarly to the voltage and result in a wave traveling through the spatial domain, the presented approach is not optimal due to the solely nonlinear model equations for the ion state probabilities and the additional costs of projecting between the reduced and full dimensional space.Instead, an approximation of the nonlinearity, for example, by the DEIM [11] can be adapted to the setting of the dynamically transformed basis functions, similar to [7, Chapter 3.2].Furthermore, since the electrophysiological model is embedded in a multiscale skeletal muscle simulation, including the contraction of the muscle, the ROM can be integrated into this much more complex framework.

Code Availability:
The Python code for the simulations is available as supplementary material and can be obtained under DOI: 10.5281/zenodo.8282665.

Figure 1 -
Figure 1 -Illustration of the space-time cylinder of the AP wave propagation for a muscle fiber of length 12 cm, originating in the center of the fiber.FOM simulation with the parameters from Table1and h t = 0.001.

3 Figure 2 -
Figure 2 -States for a time range of [0, 22] visualized from red to green with 1 ms between the curves.The first and last states are visualized in black.

Figure 3 -
Figure 3 -Spatial parts ϕ of reduced basis functions from sPOD for completely separated wavefronts initialized with the snapshot data.

Figure 4 -
Figure 4 -Absolute error of sPOD for completely separated wavefronts with r = 2, initialized with snapshot data, and using the equilibrium transformation over the short time interval T 0 , T 0 + 4 .

Figure 5 -
Figure5 -Spatial parts ϕ of reduced basis functions from sPOD for completely separated wavefronts with r = 2, zeroinitialized for the undivided and divided domain case.
captured accurately.However, for a divided domain in Figure5(b), the shape of the wavefront is reconstructed properly.

Figure 6 -
Figure 6 -Absolute error of sPOD for partly overlapping wavefronts with r = 2, initialized with snapshot data, over the short time interval T 0 , T 0 + 4 .

Figure 7 -
Figure 7 -Comparison of the original wavefront (solid line) and reconstructed wavefront (dots) of sPOD for partly overlapping wavefronts with r = 2, initialized with snapshot data for the short time interval T 0 , T 0 + 4 .

Figure 8 -
Figure 8 -Spatial parts ϕ of reduced basis functions from sPOD for partly overlapping wavefronts with r = 2 with zero initialization for the divided parts separately.

Figure 9 -
Figure 9 -Absolute error of the ROM for completely separated wavefronts simulated with h t = 0.001.

Figure 10 -
Figure 10 -Comparison of the original wavefront (solid line) and reconstructed wavefront (dots) of the ROM for completely separated wavefronts simulated with h t = 0.001 over the whole time interval [T 0 , T ].
t = 9.25 ms t = 13.5 ms t = 17.75 ms t = 22.0 ms (a) Zoom in left peaks

Figure 11 -
Figure 11 -Comparison of the original wavefront (solid line) and reconstructed wavefront (dots) of ROM for partly overlapping wavefronts simulated with h t = 0.001 over the whole time interval [T 0 , T ].

y 3 Figure 12 -
Figure 12 -Absolute error snapshots for the ion channel probability states of the ROM for completely separated wavefronts simulated with h t = 0.001.
deke (IANS and SC Simtech, University of Stuttgart) and Dr. Benjamin Unger (SC SimTech, University of Stuttgart) for the supervision of the thesis work and their continued advice.This work has been funded in parts by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC 2075/1 -390740016

Table 1 -
Parameters of the electrophysiological simulations ℓ [cm] h s [cm] T 0 [ms] T [ms] C m µF/cm 2 A = 0.001.To activate the propagation of the AP we initialize the voltage state by a scaled and shifted Gaussian function m [1/cm] σ eff [mS/cm] V eql [mV] P Na + P K +

Table 2 -
Errors of sPOD for completely separated wavefronts with r = 2, with equilibrium shift for snapshot initialized and zero-initialized for the undivided and divided domain case

Table 3 -
Errors of sPOD for partly overlapping wavefronts with r = 2 for different initializations

Table 4 -
Offline errors of the reduced ansatz space with r = 2 and runtime of the full order simulation

Table 5 -
Online errors and performance analysis of the ROM for completely separated wavefronts and various time-step widths

Table 6 -
Online errors and performance analysis of the ROM for partly overlapping wavefronts and various time-step widths