Solving a one-dimensional moving boundary problem based on wave digital principles

We report on a novel method for solving one-dimensional moving boundary problems based on wave digital principles. Here, we exploit multidimensional wave digital algorithms to derive an efﬁcient and robust algorithm for the solution of the considered problem. Our method lets the wave digital model, on which the algorithm is based, expand according to the size of the solution domain. The expanding model introduces new dynamical elements, which must be properly initialized to obtain a calculable algorithm. To deal with this problem, we make use of linear multistep methods to extrapolate the missing values. Our results show the proposed method to indeed be capable of correctly solving a one-dimensional partial differential equation describing a growing biological axon.


Introduction
Partial differential equations (PDEs) are the primary mathematical tool for describing spatiotemporal physical phenomenon, such as electro-magnetic wave propagation (Lindell, 2005) or heat transfer (Vaidya et al., 1913).In general, most partial differential equations can not be solved analytically, which is why numerical solvers are necessary.
An interesting subset of differential problems is the set of moving boundary problems also known as Stefan problems (Stefan, 1889).These problems represent expanding physical domains, which, from a mathematical perspective, implies that the solution domain is dynamical, i.e. its boundaries result as an integral part of the solution (Zerroukat & Chatwin, 1997).Such problems occur in many research fields.Examples are the analysis of melting and solidification problems (Tao, 1967;Alexiades, 1993) or diffusional controlled release problems associated with drug release (Higuchi, 1963;Lee, 2011).The solution of moving boundary problems is quite challenging, even from a numerical standpoint.Over the years, a lot of work has been invested into the development of approximative, numerical, and semi-analytical solution methods of which some turned out to be very accurate compared to real-world measurements, refer to (Zerroukat & Chatwin, 1997;Tao, 1986) for a good exposition.
In this paper, we aim to introduce a new numerical solution method based on the wave digital concept (Fettweis, 1986).To this end, we have picked one example, where such as problem occurs, namely the diffusion of microtubules in a spatially extending axon (McLean & Graham, 2004).What makes this process interesting is the fact that it plays an important role in the maturity of neural networks.From an engineering perspective, both the mathematical and electrical modeling of this process can be used for the design of self-organizing bioinspired neural networks in the context of neuromorphic computing, cf.Ochs et al. (2021), Michaelis et al. (2022), Jenderny andOchs (2022), Singh Muralidhar et al. (2022), which has inspired us to pick this specific problem.
The synthesis of a reference circuit constitutes the hardest part of the solution process.In the early stages of the wave digital concept, these circuits were derived based on a combination of physical considerations and circuit-theoretical intuition (Fettweis, 1992(Fettweis, , 2002;;Erbar & Horneber, 1995;Fettweis, 1994).However, this required a strong foundation and extensive experience in the field of circuit theory.Some circuit synthesis approaches were later presented for example in Vollmer (2004), Vollmer (2005).Later on these approaches were used as a basis for deriving an algorithm for the automated circuit synthesis of passive hyperbolic partial differential equations based on a special state-space representation (Leuer & Ochs, 2009, 2012).This paper applies the latter algorithm, while also combining numerical integration methods with the wave digital concept to solve the considered biological problem.Contrary to common approaches in literature, which utilize nonlinear time-varying transformations to deal with the moving boundary, our approach allows for retaining the simplicity of the baseline mathematical model.
The remainder of this paper is organized as follows.Section 2 recapitulates the considered problem from a biological as well as a mathematical standpoint.In Sect.3, we derive an electrical model for the considered problem.Section 4 discusses the concept of dynamic wave digital models, which is then utilized to obtain the solutions presented in Sect. 5. Finally, Sect.6 summarizes this work and its main contributions.

Biological background
Neurons or nerve cells are the primary building blocks of the human nervous systems.In general, these cells have a wide range of responsibilities including sensing environmental signals, processing sensory signals, and communicating information to different parts of the body (Bullock et al., 2005).The neurons themselves are composed of many different parts, the soma representing the cell's main body, the nucleus containing genetic information, the axon responsible for bridging the spatial gap between neurons, the synapse handling interneuron communication, etc.A simplified biological sketch of a neuron is presented in Fig. 1.In this work, we focus on one part of the neuron, namely the axon.In particular, we deal with the mathematical and electrical modeling of axonal outgrowth.
During the development of the nervous system, axon-growth takes place.The axon's form is maintained by a cytoskeleton consisting of tubulin dimers (microtubules) (McLean & Graham, 2004).Its growth process involves the extension of the cytoskeleton and the spatial elongation.Here, tubulin is first produced within the cell body and then transported along the axon (Wang et al., 1996).The latter is driven by a combination of diffusion and active transport (Galbraith et al., 1999).The axon's extension results from the polymerization of tubulin at the tip of the axon (Morrison et al., 2002), where microtubules are first assembled and then used to create the cytoskeleton, which results in axon-growth (Sayas et al., 2002).Considering that this process takes place on a timescale of months to years, we expect some form of tubulin degradation.Indeed, this effect was experimentally observed, where the degradation of tubulin took place on a timescale of days and thus slowed down the growth process (Miller & Samuels, 1997).Overall, the process of axon-growth is determined by the combination of tubulin production (at the soma), tubulin degradation, tubulin transport (along the axon via diffusion and active transport), and microtubule assembly (at the axon tip) (Miller & Samuels, 1997).
Fig. 2 A mathematical view on tubulin-driven axon-growth for three different time instants with t 0 < t 1 < t 2 .The three-dimensional cylinder approximates the axon's geometry.We assume a steady flux of tubulin being produced at the soma and flowing through the soma-axon interface into the axon.At the tip of the axon, whose current position is denoted by l(t 1 ), microtubules are assembled to form the elongation of the axon.Under the assumption of a homogenous flow of tubulin, we denote the concentration of tubulin at the coordinates (z 1 , t 1 ) by c(z 1 , t 1 )

Mathematical modeling
A mathematical view on the biological scenario of axon-growth is depicted in Fig. 2. Here, a three-dimensional cylinder is used to approximate the axon's geometry.Moreover, we assume that the axon only grows straight and in the z-direction.Here, the variable (t) denotes the axon's length at the time t.At the soma-axon interface, see Fig. 2, we assume a steady supply of tubulin leading to a constant in-flux.At the axon tip, see (t 1 ) in Fig. 2, microtubules are assembled to form the extension of the axon, so it can elongate.The axon's extension is represented by the cylinder with the black dashed edges.Under the assumption that the concentration of tubulin is identical for all coordinates (x, y) at the cross-section at z 1 , we refer to the concentration of tubulin at some time instant t 1 as c(z 1 , t 1 ).
The following section recapitulates the derivation of a one-dimensional PDE model for the tubulin distribution within the cylinder based on conservation laws.This derivation originates from the discussions of McLean and Graham (2004) and is first made for the three-dimensional case.Reducing the three-dimensional model to one dimension follows in a straightforward manner.

Conservation law
Let V a ∈ R 3 and S a ∈ R 2 denote the volume and enclosing surface of the three-dimensional cylinder of Fig. 2, respectively.Furthermore, let r = [x, y, z] T be the three-dimensional space vector.We start our derivation from the three-dimensional case, where the tubulin concentration at some coordinate r ∈ R 3 is denoted by c(r, t).Our first goal is to derive an integral conservation law from which a partial differential equation for the tubulin distribution can then be inferred.To this end, we assume the axon to be an isolated system and exploit the principle of conservation of charge from physics.Here, we can draw an analogy between the concentration of tubulin and charge density.The principle of conservation of charge states that the change in electric charge in any volume is equal to the amount of in-flowing charge minus the amount of out-flowing charge.In our context, we can similarly state that the change of tubulin within V a is equal to the flux of tubulin in and out of the axon: (1) Here, D t denotes a partial derivative w.r.t.time, j(r, t) denotes the flux of tubulin at the coordinate (r, t), and s a is the curve parametrization of S a .Equation ( 1) is an ideal conservation law, as it does not consider the possibility of tubulin degradation within the axon.We can incorporate this non-ideality by adding a degradation term, where c l (r, t) is a (possibly nonlinear) function of c(r, t) describing the loss of tubulin within V a , while T l represents the exponential decay constant, which is identical for every point r ∈ V a .Under the assumption of a smooth vector field, we can apply the divergence theorem to rewrite the conservation law (2) as where D r = [D x , D y , D z ] T denotes the gradient w.r.t. the coordinate vector r.Equation ( 3) is the conservation law in integral form.In differential form, it reads: Thus, we obtain a partial differential equation (PDE), whose solution provides the tubulin distribution in V a .However, this PDE is incomplete, since we have yet to define the tubulin flux j(r, t) and the tubulin loss function c l (r, t) in terms of c(r, t).

Tubulin flux and degradation
Within the axon, the flux of tubulin is either due to diffusion j d (r, t) or active transport j a (r, t), where the latter is meditated by transport proteins.The overall flux is given by the superposition of both processes: j(r, t) = j d (r, t) + j a (r, t) . (5) Fick's law of diffusion is arguably the most common way of describing diffusion processes (Paul et al., 2014).It states that the flux goes from the region of high concentration to the region of low concentration.Moreover, the flux is said to be proportional to the spatial derivative of the concentration: where D(r) ∈ R 3×3 is the so-called diffusivity matrix, which has a spatial dependency.
The active transport of tubulin j a (r, t) is assumed to be proportional to the local tubulin concentration: where v a (r) denotes the local transport velocity of tubulin.Inserting both definitions ( 6) and ( 7) into (5) yields Thus, we have now defined the tubulin flux j(r, t) in terms of the tubulin concentration c(r, t).
Regarding the tubulin degradation, we adopt the tubulin loss model from McLean and Graham (2004), which states that tubulin degradation is proportional to its local concentration, i.e. c l (r, t) = c(r, t).If we briefly assume a closed system (no tubulin entering or exiting the axon), then the considered form of tubulin loss implies that tubulin locally decays to 1/e of its initial value after an elapsed time of T l , i.e. c(z, t 0 ) ≈ e −1 c(z, t 0 + T l ), for some initial time t 0 .

One-dimensional model
Since we are assuming the axon to only extend in one direction, it is sufficient to work with a one-dimensional model.Moreover, both the diffusion constant and the active transport velocity are assumed to be constant along the axon, which leads to the special case: This leads to the one-dimensional PDE system where the last equation (10e) is the growth law describing the elongation or shortening of the axon.It states that the axon, with the length (t), grows or shrinks at the velocity v g , when the tubulin concentration at the axon tip c( (t), t) exceeds or drops below some constant threshold c th .For the sake of completeness, we state that 0 denotes the initial length of the axon, while N 0 = 1 mol is a normalization constant arising from the assumption that the tubulin concentration is measured in mol.
Looking at the PDE system (10), we see that the tubulin distrubiton plays a role in solving the differential equation describing the growth process (10e), but that the axon length (t) resulting from the solution of the latter equation does not directly appear in (10a)-(10d) describing the tubulin distribution.Here, it is important to consider that the solution domain of the PDE is determined by the length of the axon (t).In other words, the PDE is only solved on the domain z ∈ [0; (t)] at every time t.This makes the problem a so-called moving boundary problem, which essentially denotes a differential problem with a timevarying solution domain (Zerroukat & Chatwin, 1997).To our utmost knowledge, no general analytical solution exists for this equation.However, for some special cases in which the appearing parameters are restricted to certain intervals, McLean and Graham were able to derive analytical steady-state solutions, cf.McLean and Graham (2004).Furthermore, the existence and uniqueness of one solution for the biologically relevant case, i.e. (∞) > 0 was proven in the same reference.Moreover, we would like to emphasize that the PDE system is not hyperbolic and does not (generally) constitute a passive PDE system making it incompatible with the wave digital method of solving PDEs.We will show how the PDE system can be made hyperbolic without changing the steady-state solution.We will also derive necessary and sufficient passivity conditions.

Boundary conditions
Naturally, the considered PDE system can only be solved when boundary conditions are imposed.Here, we adopt the boundary conditions proposed by McLean and Graham (2004), which we now briefly recapitulate: 1.There is a steady (independent) flux of tubulin entering the axon at z = z 0 (Neumann-type condition).2. The flux of tubulin exiting the axon at z = z e is proportional to the tubulin concentration at the tip of the axon (Robin-type condition).
Mathematically, they can be formulated as, where c 0 is a typical tubulin scale, r p is a typical tubulin production rate, r a is the tubulin assembly rate, and j d,r is the constant returned flux due to tubulin disassembly.In biology, the latter quantity is generally not constant.To our utmost knowledge, its value stems from a steady-state measurement and is used as an approximation (McLean & Graham, 2004;Graham et al., 2006).Note, both assumptions at the boundaries describe diffusion processes, which is why the boundary conditions are formulated in terms of and not j(z, t).If we formulate the conditions in terms of j(z, t), we would also impose restrictions on the active transport of tubulin.However, the latter transport mechanism remains unchanged at the boundaries of the considered domain.

Initial conditions
In addition to the boundary conditions, we require initial values for the starting time t 0 .In the following, we use the linear initial conditions suggested in Graham et al. (2006): These initial conditions fulfill the boundary conditions (11) but not the governing PDE (10).
The initial conditions for j(z, t) can be obtained by using (10b), which, in compact form, reads: The spatial derivative of c(z, t) at t = t 0 can be obtained from (12a):

Electrical modeling
Solving PDEs with the wave digital concept requires first synthesizing a reference circuit, whose dynamics coincides with the considered PDE (Fettweis & Nitsche, 1991a).In the following, we exploit the methods discussed in Leuer and Ochs (2009), Leuer and Ochs (2012) to synthesize a reference circuit.

Physically-motivated hyperbolization
The wave digital concept is not equally suited to solve all types of PDEs (Fettweis & Nitsche, 1991a;Luhmann & Ochs, 2006).It excels at solving hyperbolic PDEs, because all other types of PDEs lead to wave digital algorithms with implicit relationships, which must be resolved, for example, through fixed-point iterations, cf.Luhmann and Ochs (2006).The implicitness of these relationships stems from the lack of deep physical modeling allowing physical phenomenon to propagate at infinite speeds.The considered PDE ( 10) is of the parabolic type.For a deeper physical modeling (and an enhanced wave digital algorithm later on), we now apply a hyperbolization by finding the root cause, which makes the PDE parabolic, and rectifying it such that the PDE becomes hyperbolic.The root cause is the use of Fick's law, which does not consider any sort of delay in the diffusion of tubulin.This implies that even for very small times, there exists a finite amount of diffusing tubulin, whose origin is very far from its current position, cf.Gómez et al. (2007).This makes the diffusion process unphysical and distorts the true qualitative nature of the considered phenomenon, as the movement of tubulin particles is exceeding the speed of light.Thus, for a deeper physical modeling, we will now replace Fick's diffusion law (10c) by a modified diffusion law, namely Cattaneo's law Cataneo (1958): Cattaneo's law is the most the widely accepted generalization of Fick's diffusion law (Gomez et al. 2010).It introduces a sort of inertia or relaxation time with the time constant T r into the diffusion process allowing the particles to only move at an infinite propagation velocity, whose value will be derived later on.Considering that our PDE system represents a general advection-diffusion problem on the most abstract level, making use of Cattaneo's law is a very common and natural way of dealing with the unphysical nature of this problem.A general theory on hyperbolic advection-diffusion problems utilizing Cattaneo's law has been developed in Gomez et al. (2010), and we refer the interested reader to this reference for more details on this hyperbolization approach and its physical consequences.Note, the underlying conservation law governing our PDE system is retained, when modifying the diffusion law, cf.Gómez et al. (2007), Gomez et al. (2010).Also, we would like to emphasize that there is a great amount of literature concerning the generalized hyperbolization of PDEs modeling physical processes besides Cattaneo's method (Fettweis & Basu, 2015;Bright & Zhang, 2009;Rubin, 1992;Chester, 1963), some more recent and/or physical than others.With reference to (10), the hyperbolized one-dimensional PDE now reads: One important property of this hyperbolic PDE system is the fact that it degenerates to the parabolic model in the steady-state, i.e. when all time derivatives amount to zero.This means that its steady-state solution does not differ from the one of the parabolic model, such that the existence and uniqueness properties of its solution (McLean & Graham, 2004) carries over from the parabolic model.For the remainder of this section, we consider the two parts of (14) as two independent systems.This is because the c(z, t) appears in the ODE computing (t), but (t) does not appear in the PDE computing c(z, t).Thus, we consider the ODE to be controlled by the PDE and therefore first start by synthesizing a reference circuit for the PDE before dealing with the ODE.

Canonical PDE formulation
Before applying the circuit synthesis procedure, one must first rewrite (14a) in the canonical form, where the voltage vector e is a vector of zeros due to the absence of source terms in the left part of ( 14).The relation between the current vector i and the tubulin concentration and flux is given by: Here, we introduce the normalization constants F 0 = 1/s, L 0 = 1 H, and associate c(z, t) and j(z, t) with currents.Now, the matrices L t , L z , and R are obtained by first rewriting ( 14) in terms of the currents defined in (15b) and then applying a coefficient comparison: Note, the circuit synthesis procedure only works when the PDE fulfills the conditions of multidimensional passivity (MD passivity) (Fettweis & Basu, 2011;Leuer & Ochs, 2009, 2012).The latter is given when the circuit is MD causal (Fettweis & Basu, 2011) and the following sufficient conditions are fulfilled by the involved matrices: We elaborate on the conditions of MD causality in Sect.3.5.2.Now, we have decomposed the PDE so the first two conditions are met.The third condition is generally not fulfilled, since the matrix R is only positive-semidefinite for: This means that the circuit synthesis can result in an active (resistive) component.For the sake of completeness, we would like to mention that the parameters used in this manuscript violate this condition, which then indeed results in an active component.In this case, many of the great properties inherent to multidimensional wave digital algorithms may or may not hold.We will come back to this aspect in Sect.3.6, once we complete the circuit synthesis procedure.

Circuit realization of the axon-growth PDE
The canonical formulation (15) can now be used as a starting point the systematic circuit synthesis.In the following, we only briefly sketch the procedure; the interested reader is referred to Leuer and Ochs (2009), Leuer and Ochs (2012).To this end, we first rewrite (15a) as: Here, the left and right handside are termed (hyperbolic) PDE and ODE subsystem, respectively: To ensure the hyperbolicity of the PDE subsystem, a new term L z t D t has been introduced.The matrix L z t is a positive-semidefinite diagonal matrix consisting of free/design parameters with the unit of an inductance: The matrix is defined such that a positive free parameter α z μ is introduced every time the matrix L z has a non-zero column (or row, since it is symmetric).This ensures that a time derivative is added to every spatial derivative introduced by the term L z D z , which guarantees the hyperbolicity of the PDE subsystem.When no spatial derivative appears in a certain row of the equation set, i.e. when the corresponding column fulfills L z e μ = 0, no time derivative is required and the corresponding free parameter α z μ is set to zero.Due to the introduction of L z t and because the original PDE (15) should not be modified, one must introduce a new matrix L t t , which is defined as: By inserting (18e) into (18a), it can be verified that we indeed re-obtain the canonical formulation (15).The newly introduced matrices are given by: For the sake of completeness, we now give a brief outline on the circuit synthesis procedure proposed by Leuer and Ochs (2009), Leuer and Ochs (2012).First, to ensure MD causality, the hyperbolic PDE subsystem (18b) undergoes a bijective coordinate transformation, cf.Fettweis and Basu (2011).Second, the transformed PDE subsystem is realized with a known equivalent circuit.This results in an uncoupled circuit representing the (uncoupled) PDE subsystem.Third, an equivalent circuit for the ODE subsystem is synthesized by applying the circuit synthesis procedure extensively described in Ochs (2012).Fourth, the circuit realization of the ODE subsystem is used to close the external ports of the circuit representing the PDE subsystem.This step is a direct consequence of the constitutive relationship (18a).Finally, the previously introduced free parameters α z μν , can be chosen in a way to eliminate redundant electrical components and obtain a simpler circuit realization.This step can be understood as a circuit optimization.Applying all aforementioned steps yields the reference circuit C 1 in Fig. 3 with the electrical parameters, where the newly introduced velocity constant v 1 is used as a transformation variable, cf.Fettweis and Nitsche (1991a), Leuer and Ochs (2012).In many cases, an appropriate choice of v 1 can greatly simplify the resulting circuit.We elaborate on the optimal choice of this constant in Sect.3.5.Note, we refer to the inductors with the derivative operators D z ± as spatio-temporal inductors, because the derivative operator is given by a linear combination of a spatial and a time derivative, cf.Fettweis and Nitsche (1991a), Luhmann and Ochs (2006), Leuer and Ochs (2009), Leuer and Ochs (2012): As can be seen from the definition of R 1 in ( 19), its resistance is only non-negative when ( 17) is fulfilled.In other words, the active part of the PDE translates into a negative resistance in the resulting circuit, which does not dissipate but rather supplies power to the circuit.
While the applied procedure indeed returns a valid reference circuit, the number of components is quite large in comparison to the order of the PDE.For that reason, we now adopt a different approach for obtaining a reference circuit.To be specific, we use the systematic procedure only to synthesize a reference circuit for the PDE subsystem.By doing so, we obtain the reference circuit depicted in the dashed box Fig. 3. To obtain a reference circuit for the remaining ODE subsystem, we first explicitly write out the associated equation (18c): Sorting the terms and neglecting the appearing zero vector, we obtain the two mesh equations: Let us first consider the first mesh equation.We see that it corresponds to the series connection of an inductor L 1 , a resistor r 1 , and the port with the voltage u 1 in C 2 of Fig. 3. Now, let us consider the second mesh equation, while ignoring its right handside.This equation is identical to first mesh equation, and we effectively obtain the same electrical components with different parameters, which are interconnected to the port designated with the voltage u 2 in Fig. 3. Contrary to the first equation, the second equation contains a source term r 21 i 1 coupling i 1 to the other side of the circuit.To realize the cross-coupling, we can make use of a current-controlled voltage source with the transresistance r 21 encoding the active part of the PDE.The reference circuit resulting from this approach is depicted in C 2 of Fig. 3.As can be seen, it contains less electrical components and has a simpler structure than C 1 in Fig. 3.
For the remainder of this work, we make use of the alternative circuit, due to its structural simplicity.

Circuit realization of the axon-growth ODE
In this section, we derive an equivalent circuit for the growth dynamics given in (14b).To this end, we associate the axon's length with a current quantity: Then, we substitute the definition of D t (t) from (10e) rewrite its right handside in terms of the currents defined in (15b): As can be seen, ( 22b) performs an integration of its right handside, when (t) is known.Thus, a possible circuit realization is given by the simple integrator circuit depicted in C 4 of Fig. 3.
Here, r t denotes the transresistance of the appearing current-controlled voltage source.The axon's length can be reconstructed from the current i 3 flowing through the inductor L a by applying the bijective mapping relation given in (22a).

Circuit optimization
The last step of the circuit synthesis procedure involves optimizing the circuit C 2 in Fig. 3 by eliminating redundant components, which can be achieved, by an appropriate choice of the free parameters α z μν and the velocity constant v 1 .Here, we start by determining the maximal propagation velocity v max of the system, which allows us to derive a lower bound for the velocity constant v 1 (Fettweis & Nitsche, 1991a;Gomez et al., 2010).

Maximal propagation velocity
For hyperbolic PDEs with constant coefficients of the form (18), the maximal propagation velocity v max can be explicitly calculated by evaluating the dispersion relation (Bilbao, 2004), where det(•) is the determinant operator, j is the imaginary unit, and k z is the wavenumber in z-direction.The dispersion relation provides the angular frequencies ω ν of the waves, which propagate in the system, from which we can derive the associated group velocities: The maximal propagation velocity is the maximal (absolute) group velocity, which turns out to be a function of k z and not a constant.A variable maximal propagation velocity is not helpful for the upcoming circuit optimization procedure, which is why we now make a worst case approximation.Considering that (24a) is a monotonically increasing uneven bounded function of k z , we can make a worst case approximation by calculating the limit:

Simplifications
To eliminate redundancies in the circuit, we now choose our free parameters such that certain electrical parameters become zero.Considering that our two free parameters can be chosen to eliminate the inductances in ( 19), we conclude that the following two choices are optimal: The two design parameters have influence on the other electrical parameters of the PDE subsystem: The final parametric degree of freedom is given by the velocity constant v 1 .This constant must fulfill the MD causality condition (Fettweis & Basu, 2011;Leuer & Ochs, 2009), which in this case reads: v 1 ≥ v max .Furthermore, it must be chosen so the inductances L z 12 and L z 21 are non-negative, which is a necessary condition for their passivity: Since v max > 0, the second condition is actually redundant.The non-negativity condition of L z 12 coincides with the condition of MD causality.Thus, we see that v 1 = v max is an optimal choice, since it eliminates L z 12 without violating the latter.In a final optimization step, we eliminate the ideal transformer, appearing in C 2 of Fig. 3, by an equivalency transformation incorporating its turns ratio into the electrical parameters on the left side of its primary port.The reference circuit resulting after these simplifications is depicted in C 3 of Fig. 3.The equivalency of its dynamics with those of PDE system can be verified by evaluating the node and mesh equations.

Intermediate summary and discussion
In an intermediate summary, we would like to briefly elaborate on the properties of the synthesized circuit and its consequences on the resulting wave digital algorithm.In this work, we have synthesized two possible electrical circuits for the hyperbolic PDE system.If the circuit parameters are chosen in a way that fulfills the third passivity condition (17), we recommend using C 1 as depicted at the top of Fig. 3 to solve the hyperbolic PDE.However, in case this condition is not met, we recommend using circuit C 3 depicted at the bottom of the same figure due to its structural simplicity.The latter circuit contains a nonlinear and active component and as mentioned earlier, there is no guarantee that any of the great properties inherent to multidimensional wave digital algorithm carry over to this case.However, we can see that all circuit elements are passive and that the active part of the PDE is reflected by one active nonlinearity, the controlled voltage source.Thus, we can separate the circuit into a passive subcircuit that is interconnected to an active subcircuit, a controlled voltage source.All great properties of wave digital algorithms apply to the passive and dominant part of the overall circuit.In this configuration, the active subcircuit excites the passive subcircuit.As long as this excitation is bounded the passive subcircuit will react with a bounded output due to passivity.Here, it is important to consider that passivity carries over to the wave digital domain (Meerkötter, 2018).

Wave digital modeling
The emulation of the electrical model derived in the previous section is discussed in this section.For a better understanding of the applied methods, we begin by briefly recapitulating the fundamentals of the wave digital concept.Then, we translate the electrical model into the wave digital domain, which results in a compact wave digital model.Uncoiling the latter results in the explicit wave digital model, which can be used for emulation purposes.Finally, we introduce the concept of dynamic wave digital models, which is our main contribution.

Fundamentals
The wave digital algorithm of a given reference circuit is the result of writing out the signal relationships of the associated wave digital model in a causal manner (Fettweis, 1986;Meerkötter, 2018).The latter can be obtained by decomposing the reference circuit into a set of one-and multiports, which are then translated into the wave digital domain by applying the bijective mapping relation: Here, a, b, and R denote the incident wave, reflected wave, and the port resistance, respectively (Fettweis, 1986).Electrical components with differential relationships must first undergo a discretization via numerical integration, commonly the trapezoidal rule (Fettweis, 1986;Meerkötter, 2018).For an overview on electrical components and their translation into the wave digital domain, the interested reader is referred to Fettweis (1986), Meerkötter (2018).
Considering that we will be working with multidimensional wave digital algorithms, we define the discrete coordinates, z k = z 0 + k Z, where z 0 denotes the initial z-coordinate.Thus, we have z k − z k−1 = Z for all k ∈ N + 0 , where T and Z = v 1 T denote the temporal and spatial step-size of the wave digital algorithm, respectively.
One aspect we would like to highlight is the fact that the trapezoidal rule, used in combination with the wave digital concept, is a so-called A-stable integration method (Hairer & Wanner, 1991).This means that asymptotic stability properties exhibited by the solution of a differential equation carry over to the wave digital domain.Indeed, this special property of the trapezoidal rule has been shown to play a big role in the accuracy of numerical solutions, for example, more recently in the simulation of lumped transmission lines, cf.Zhou et al. (2021).In our context, this is property also plays a big role, as the existence and uniqueness of a steady-state solution is ensured by the discussions of McLean and Graham (2004).Thus, multidimensional wave digital algorithms based on the trapezoidal rule will always have an edge over any explicit integration method for solving this problem for two reasons: 1.The trapezoidal rule is an implicit, lossless, and the most accurate A-stable linear multistep method (Dahlquist, 1963;Ochs, 2001), which becomes explicit in the wave digital domain (Fettweis, 1986).Consequently, the wave digital algorithm is fully explicit and does not require any sort of iteration methods, which makes it efficient.2. Explicit integration methods can not be A-stable (Dahlquist, 1963;Ochs, 2001).Contrary to A-stable integration methods, explicit integration methods are only numerically stable for sufficiently small step-sizes (Dahlquist, 1963), which can make them inefficient from a computational perspective.

Wave digital model
Now, we discuss the translation of the circuits C 3/4 in Fig. 3 into the wave digital domain.Starting with circuit C 4 , the ideal controlled voltage source and the inductor translate to a reflective wave source with the reflection coefficient −1 and a delay element (with a sign inversion), respectively (Fettweis, 1986).The port resistance of the temporal inductor is given by R a = 2L a /T , where T denotes the sampling period of the wave digital algorithm.Now, we translate circuit C 3 , where we start by replacing the Jaumann structure with a Jaumann adaptor.The resistor and the controlled voltage source are then replaced by the reflection coefficient R and a reflective wave source, respectively, with: Finally, the spatio-temporal inductors are replaced by spatio-temporal delays Z ± (Leuer & Ochs, 2009, 2012).The delays are associated with transformed differential operators ( 20) and can be interpreted as a combination of a time and a positive or negative spatial delay (Luhmann & Ochs, 2006).To understand the reason behind two signs of the spatial delay, one must consider that physical phenomenon can propagate in positive or negative directions as time passes.
The compact wave digital model M 1 is depicted in Fig. 4. It is an implicit representation of the actual wave digital model, because it only describes the dynamics of the PDE in one grid point.The uncoiled wave digital model, on the other hand, is an explicit interpretation of a multidimensional wave digital model (Luhmann & Ochs, 2006), see M 3 in Fig. 4. Since we are working with a one-dimensional model, consider a one-dimensional grid with spatial step-size Z .Now, place the compact wave digital model of the previous section at every grid point and remove all the spatio-temporal delays Z ± .The ports at which the latter elements were previously placed, are left open (as external docking ports).The uncoiled wave digital model is obtained by interconnecting the compact wave digital models at every grid point.In particular, every reactive-element-free model is only connected to its nearest neighbor.The interconnection elements are the delays, which we have previously removed from the model.These delays, however, are now pure temporal delays.This is because the spatial delay of the compact model is now embedded into the distribution of the compact models along the one-dimensional grid.For instance, the currents i ν (z 1 , t) can be inferred from the wave quantities of the first wave digital model (the one placed at the first grid point) in M 3 of Fig. 4. When interconnecting two neighboring models, it is important that every connection exiting an external port, which previously belonged to a positive/negative spatio-temporal delay, also connects to a port belonging to the same type of delay.Note, in this work, we do not discuss the realization of the boundary conditions.However, our approach is based on the discussions of Fettweis (2006).

Dynamic wave digital model
In the first section, we stated that the considered problem is of the moving boundary type.In this section, we first start by explaining the consequences of a moving boundary w.r.t.wave digital model and then move on to demonstrate a way of dealing with these consequences.When an axon changes its length, the tip of the axon, which was previously at some position (t 0 ) = 0 , moves to some position (t 1 ) = 1 at the next time instant, where 0 can be larger than, smaller than, or equal to 1 .

Fig. 4
Compact representation of the PDE subsystem in the wave digital domain (M 1 ) C 3 .The wave digital model M 2 corresponds to translating the equivalent circuit of the growth law (10e), denoted with C 4 in Fig. 3, to the wave digital domain.Uncoiling the compact wave digital model M 1 leads to the uncoiled wave digital model M 3 , in the case of three spatial coordinates z μ with μ ∈ {1, 2, 3} Case 1: Growing axon First, let us consider the case, where 0 is larger than 1 .At the time step t 0 , the right boundary of the solution domain is at the position 0 .In the next time step, the axon grows, such that the boundary is now at 1 .Thus, the position at which the (right) boundary condition must be applied, has changed.In terms of the wave digital model, this means that we must add a wave digital structure to the wave digital model every time the axon grows.Furthermore, we must shift the position of the right boundary to be at the new boundary coordinate, which comes after the newly introduced wave digital structure.Here, we have implicitly assumed the difference between the new axon length and the old axon length to be maximally given by Z .Every time a wave digital structure is added, two new delays are introduced, see the right side of Fig. 5.The delays couple the wave digital model from the previous time step to the newly added wave digital structure.To initialize the delay that the models have in common (at the two consecutive time steps), we can simply use the wave reflected towards the right boundary at the previous time step.However, we have no way of initializing the other introduced delays, such that, our wave digital model becomes incalculable.To deal with this problem, we suggest using an extrapolation technique to calculate the missing values, which is the topic of the next section.Case 2: Shrinking axon Now that we have covered the first case, let us discuss the second case, where the axon shrinks, i.e. 1 < 0 .In this case, we can simply wipe the memory of all the wave digital structures representing the axon segments that come after 1 .Now, we have a wave digital model that is fit to solve the PDE in the domain [z 0 ; 1 ].Essentially, this corresponds to interpreting Fig. 5 in the reverse manner, i.e. the right side is the initial wave digital model and the left side is the wave digital model at the next time instant.

Case 3: Axon retains length
The third and final case is when the axon retains its length.This case is very trivial, as we can simply continue using the same wave digital model to solve the PDE on the same solution domain, until one of the other two cases comes into play.

Extrapolation technique
Extrapolation techniques are useful for approximating unknown values from known values within the solution domain.In this section, we focus on a special type of extrapolation techniques, namely linear multistep methods (LMS), cf.Fränken and Ochs (2002) and derive a LMS method for approximating the initial value of the newly introduced delays.These values depend on the tubulin concentration c(z, t) and the tubulin flux j(z, t).Thus, the initialization requires extrapolating the associated currents i 1 and i 2 .In the following, we demonstrate the extrapolation for the missing current i 1 .The derived method can then be analogously used to extrapolate i 2 .
Let us consider the current i 1 (z, t) associated with the tubulin concentration c(z, t) at some point z k .We can approximate the currents at the preceding coordinates i 1 (z k−ν , t) by a Taylor series expansion, while using i 1 (z k , t) as our pivot point: Here, the notation D z c(z k , t) refers to the spatial derivative of c(z, t) evaluated at z k .Using z k−ν − z k = −ν Z and repeating this expansion for ν = {1, 2, 3} yields the linear equation Let us now first consider the left side of the above equation.The coefficient matrix on the left handside is a so-called Vandermonde-matrix.Its inverse can always be explicitly calculated.
For t = t μ , the vector on the right handside contains known current values from the solution domain.The vector of unknowns contains the boundary current i 1 (z k , t) and its scaled spatial derivatives.To calculate the value of i 1 (z k , t μ ), we can simply use the inverse of the coefficient matrix at every t μ .Now, consider the right part of ( 28).Since we are primarily interested in calculating i 1 (z k , t μ ), we only require the first row of the equation system, which reads: Equation ( 30a) is a second-order accurate approximation for the value of i 1 (z k , t μ ).This formula also applies when extrapolating the current i 2 (z k , t μ ) associated with the tubulin flux.Using these two values, we can approximate the missing currents, which we now associate with the new wave digital structure preceding the boundary, see right side of Fig. 5. Using the half-step approximation scheme (Hetmanczyk & Ochs, 2009, 2011) and the extrapolated currents, we are now able to initialize the new delays.We can approximate the internal states of the delay elements by knowing the corresponding currents half a step in the past, cf.Hetmanczyk and Ochs (2009), Hetmanczyk and Ochs (2011): The extrapolated currents serve as an approximation for the unknown currents and prove to be sufficiently good for solving the considered PDE.Once the delays are initialized, the wave flow diagram must be recalculated before moving to the next time step.This means that two wave flow diagram iterations must be executed every time axon-growth takes place.This re-iteration is necessary to calculate the new value that is reflected by the right boundary, since axon-growth has caused it to shift its position to the right.
At first glance, it may be appealing to use a higher order approximation to calculate the missing values at the boundary.In fact, calculating the inverse of the Vandermonde-matrix would not be difficult even if we used all n − 2 points within the solution domain.Moreover, the inverse would only have to be calculated once at the beginning of the simulation, so it would not lead to any decrease in the algorithmic efficiency.However, one must not forget that the accuracy of the Taylor series expansion decreases the further we move from the pivot point.These inaccuracies are so large such that even fourth-order approximations can lead to instable solutions.Thus, we recommend the use of second-order or third-order approximation schemes, when it comes to dynamic model expansion.
Finally, we would like to mention that our method only works well under the assumption made in Sect.4.3, namely, that the difference between boundary position at two consecutive instants is maximally given by Z .In many physical problems, this assumption is justified and can be ensured by choosing the sampling period T in a way that allows the algorithm to  10) and ( 11 closely capture the movement of the boundary.In case the boundary growth is rapid, such that the relative increase in the boundary position is more than Z , the only solution would be to modify the growth law such that the boundary grows at a slower rate.Otherwise, the extrapolation method would lead to large numerical errors possibly leading to inaccurate solutions. In this scenario, the length of the axon is mainly determined by the production rate r p .Here, it was reported that the axon reaches a steady-state axon length of ∞ ≈ 35 mm for the given production rate.This results in a phenomenon, which is referred to as moderate growth, i.e. the steady-state axon length is moderately larger than its initial length.In Graham et al. (2006), the steady-state length was shown to be slightly under 1 mm.Furthermore, the peak axon length was approximately 4 mm.The steady-state concentrations at the two tips of the axon were shown to have the values c s (∞) ≈ 10.2 µmol and c t (∞) ≈ 10 µmol .
Our results depicted in Fig. 6 match those from Graham et al. (2006).The peak axon length and the steady-state length are the same.However, the tubulin concentration at the soma, which in our case amounts to c s (∞) ≈ 10.16 µmol is 0.004% smaller than in Graham et al. (2006).Furthermore, the characteristic behavior of the axon length, the tubulin concentration at the soma and axon tip, and the steady-state tubulin concentration are all the same.

Sensitivity analysis
Now, we exploit the wave digital concept to perform a simple sensitivity analysis, where we (independently) vary the diffusion constant d, the active transportation velocity v a , and the tubulin loss time constant T l .This analysis is inspired by the sensitivity analysis performed in Graham et al. (2006).Its main purpose is to verify the validity of our model and compare its dynamical behavior to the parabolic model.Due to the absence of reference numerical data, we can only give a qualitative comparison between our results and the ones from literature. Figure 7 depicts the results of our sensitivity analysis.Here, we examine the axon-growth dynamics for the three different growth scenarios: large (r p = 40 1 m ), moderate (r p = 20 1 m ), and small (r p = 2 1 m ).Each column depicts the results of varying the three different parameters for one growth scenario.

Large growth
The sensitivity analysis results for the large growth scenario are depicted in the first column of Fig. 7. Comparing these results to the ones from literature, we have found our results to be nearly identical.Some slight differences can be observed in the top plot, where we vary d.In Graham et al. (2006), varying the diffusion constant has no effect on the curve characteristics or the steady-state length whatsoever.In our case, the convergence speed is effected, i.e. lowering the diffusion constant makes the axon length converge slower.This observation is reasonable, since we have replaced Fick's law with Cattaneo's law, which introduced a time derivative and can hence slow down the speed of convergence.Varying the other two parameters does not seem to lead to any deviations, just like in literature.

Moderate growth
The sensitivity analysis results for the moderate growth scenario is depicted in the middle column of Fig. 7. Comparing our results to the ones from literature, we can observe very large deviations in the transient behavior.While the curve characteristics of our results match the ones from literature, prominent points such as the peak value of the axon length are different.The greatest deviation can be observed in the upper plot, where we varied the diffusion constant d.In Graham et al. (2006), varying the diffusion constant did not change the peak value of the axon length nor the convergence speed.However, it changed the steadystate length.In our results, we can see that varying the diffusion constant changes both the convergence behavior and the peak values of the axon length.In particular, reducing d seems to lead to fast convergence and lower peak values.The steady-state axon length, however, does

Small growth
The last set of results for the small growth scenario is depicted in the right column of Fig. 7. Overall, the results are very similar to the ones from literature.However, it is noteworthy to state that the hyperbolization seems to reduce the system's oscillation tendency.In particular, the axon length does not seem to overshoot as much before converging to its steady-state value, when compared to the results in Graham et al. (2006).Discussion Overall, our sensitivity results greatly match the ones from literature.In particular, our steadystate solutions are identical, which is to be expected, as we have shown the hyperbolization to not affect them.The largest discrepancies can be observed in the first row in Fig. 7, where we have varied the diffusion constant d.Otherwise, all other results match the reference results, which have been qualitatively justified in Graham et al. (2006) by experimental data.We would like to mention that the authors had no access to these measurements, which is why a direct comparison was not possible.In the following, we discuss the discrepancies appearing in the first row of Fig. 7 from a physical and biological perspective.
The physical perspective.As stated in Sect.3.1, hyperbolizing the PDE using Cattaneo's law essentially corresponds to adding a sort of inertia, which, in our context, limits the velocity at which tubulin particles can move within the axon.The original (parabolic) model allows the tubulin particles to move at infinite velocities, which is unphysical.Looking at the results of Sect.3.5.1,we see that the relaxation constant T r , contributes to setting a maximal propagation velocity of tubulin, see the definition of v max .Thus, we argue that our model serves as a better approximation of the real tubulin diffusion process, which occurs at the maximal velocity v max .
The biological perspective.If we briefly assume the results of the parabolic model to be correct, then this implies that varying the diffusivity does not influence the peak axon length for a constant tubulin production rate.However, the diffusivity of tubulin is what effectively determines the amount of particles reaching the axon tip per unit time.Therefore, the old set of results implies that the same amount of microtubules can reach the tip of the axon in the same amount of time, even when the diffusivity is decreased.This statement is not reasonable, since a lower diffusivity should also decrease the amount of microtubules reaching the axon tip per unit time.Thus, a lower diffusivity should also lead to a lower peak axon length according to the growth law (10e), which is something that can be observed in the first row of Fig. 7. Therefore, we argue that our model is more consistent with the modeled biological phenomenon.
A question that remains to be answered is how the relaxation time constant T r , introduced by our hyperbolization approach, is to be chosen relative to the growth time constant T a .In this work, we have chosen the former as T r = 1 s, which is in the same order as T a .Decreasing the relaxation time constant leads to results that are more similar to the ones from literature, which, according to our discussions, do not fully reflect the actual biological scenario.To fit the emulation results to actual measurements, this degree of freedom must be chosen appropriately, since it is directly correlated with the maximal velocity of the diffusion process.As this work is mainly concerned with the electrical modeling and wave digital emulation of the considered PDE, all discussions related to this topic are omitted.

Conclusion
In this work, we have dealt with the solution of moving boundary problems based on wave digital principles.Here, we first started by describing the problem of axon-growth from a biological as well as a mathematical standpoint.The latter is given by a one-dimensional partial differential equation with a time-varying solution domain.Then, we derived an electrical model (reference circuit) for the partial differential equation from which we obtained the associated wave digital model.To deal with the moving boundary, we introduced the concept of dynamic wave digital models, which can red their size according to the solution domain.This concept was then used to solve the considered partial differential equation.Here, we showed an astounding resemblance between our results and those from previous research, where the mathematical model originates from.
In the future, we aim to extend this method to solve two-and three-dimensional problems.This works will serve as a basis for deriving novel methods for dealing with higher dimensional moving boundary problems.Furthermore, we aim to integrate the linear multistep methods used in this work into the wave digital concept in order to obtain a unitary wave digital method for dealing with general moving boundary problems.

Fig. 1
Fig. 1 Sketch of a biological presynaptic neuron connected to the dendrites of a postsynaptic neuron.All the main components of the presynaptic neuron are indicated

Fig. 3
Fig.3Top: Reference circuit resulting from applying the circuit synthesis procedure.Center: An alternative reference circuit equivalent to C 1 .Bottom left: Reference circuit C 2 after applying simplifications by an appropriate choice of the free parameters α μ .Bottom right: Reference circuit of the axon growth law (10e)

Fig. 5
Fig. 5Left: Representation of the wave digital model at a time instant t μ .The multiports N μ represent the compact wave digital models after removing the spatio-temporal delays.Right: Axon-growth takes place and a wave digital structure is added at the time instant t μ+1 .Adding a wave digital structure introduces three new delays.Furthermore, the right boundary moves N r and docks to the newly introduced wave digital structure at z k

Fig. 6
Fig. 6 Emulation results for the biological parameters listed in table 1. Top left: Axon length over time.Top right: tubulin concentration at soma over time.Bottom left: tubulin concentration at the axon tip over time.Bottom right: Steady-state tubulin concentration over space

Fig. 7
Fig. 7 Sensitivity analysis of the axon-growth dynamics.Each column presents the axon growth dynamics for a different tubulin production rate r p .Only one parameter is varied in every row

Table 1
Biological parameters appearing in (