Velocity–Amplitude Relationship in the Gray–Scott Autowave Model in Isolated Conditions

Velocity and amplitude are two basic characteristics of any autowave, and their relationship reflects the internal regulation of the autowave system. This study proposes an approach to approximately estimate steady velocity–amplitude (VA) relation without deriving separate formulas for V and A. The approach presumes constructing an ansatz which represents the “petal” form of phase trajectory and contains V, A, and a free parameter (parameters). After substituting this ansatz, integration of model equations leads to a VA relation analytically. A free parameter (parameters) can be determined by comparing the analytical VA relation to the numerical data. As an example, we used the simplest autowave model possessing threshold, that is, the Gray–Scott model (cubic autocatalysis with linear inhibition) in isolated conditions with an immobilized precursor and a diffusible reactant. For all values of the inhibition rate constant allowing autowave solution (i.e., except approaching zero), the free parameter of ansatz belongs to a narrow interval has little effect on VA relation and can be regarded as fixed. Assumption of precursor immobilization does not influence the results qualitatively. This approach will be useful in investigations of more complex active media systems in biochemistry, combustion, and disease control.


INTRODUCTION
Active media can be generated from coupling positive feedbacks with mass transfer. They belong to the list of the most intriguing themes in chemical kinetics 1−4 and biology. 3−7 One of the most potent ways to theoretically investigate continuous active media systems is to use reaction−diffusion (RD) mathematical models of the least possible complexity. Such models give insights into the basic, qualitative properties of the system. Low-dimensional models can be either obtained by reducing large-scale detailed ones (usually comprising dozens of variables) or constructing phenomenologically from known/assumed general properties and/or (bio)chemical organization of the system. In this study, we focus on a special class of solutions of such equations, using the term "autowave" to refer to an active wave travelling in the form of self-sustaining pulse (with front, top and back), to distinguish it from a travelling front (backless).
Two basic characteristics of an autowave are velocity and amplitude. Generally, velocity and amplitude differently depend on parameters of reactions and diffusion. Thus, different systems should have different velocity−amplitude (VA) relations, and VA relation may be considered as a characteristic of the internal regulation of the system. This view poses the question of how VA relation can be derived from model equations.
The specific form of kinetic terms in RD equations (namely, nullclines' shape, number and types of steady-states) largely determines properties of the system and suitable investigation approaches. 1−10 Nonlinearity of equations usually prevents their explicit analytical integration except for the simplest models which have travelling front solutions such as reduced Nagumo (or, purely cubic autocatalysis) and Kolmogorov− Petrovskii−Piskunov−Fisher (quadratic autocatalysis) models. 1,2,[5][6][7]10 Interestingly, in the latter case, the explicit solution is known only for the special front velocity. 7,11 Even slightly more complex models are usually studied numerically or approximately. Approximate solutions are acceptable in mathematical biology and biophysics because models and their parameters are not exact, and qualitative properties of solutions are of greater interest.
Various analytical approaches have been developed for obtaining approximate solutions of RD equations and studying their properties: narrow reaction zone method (neglecting some reaction and diffusion terms in different space regions), 1 piecewise-linear approximation (replacing non-linear kinetic terms by some piecewise-linear function), 5−7 singular perturbation, or time/length scale separation analysis (investigating the limits of small parameter(s) tending to zero), 4,5,8,9 using ansatz (guessing the structure of a solution either in space or in a phase plane) 2,5,7,10,11 and their combinations. In most cases, these approaches were applied to study travelling fronts and autowaves with saturated amplitude. So far, general methods to estimate the VA relation for autowaves do not exist, to the best of our knowledge.
In this article, we propose the new, to the best of our knowledge, analytical−numerical method for estimating the VA relation in autowave systems. The method is based on approximating the phase trajectory of the system by an ansatzsingle smooth function rendering solutions in all space regions and containing V, A, and a free parameter (parameters). After substituting this ansatz, integration of model equations leads to a VA relation analytically. A free parameter (parameters) of the ansatz can be obtained by comparing this VA relation and direct numerical solution data. As an example, we take the simplest autowave model with threshold properties: a two-component system with precursor (v) conversion to the autocatalyst (u), obeying cubic kinetics, followed by a linear inhibition of the autocatalyst, see Scheme 1. Sometimes it is referred to as the Gray−Scott (GS) model. We suppose that the system is isolated (i.e., the influx of components is absent).
Various complex systems can be studied after reducing their kinetic scheme to the GS model. 2 Examples include glycolysis, 12 catalytic chemical reactors and cold flame propagation, 1,9,13−17 fungal mycelia growth, 18 and disease spreading. 19 However, available studies in the spatialdistributed conditions are almost always limited to open conditions, 9,17−20 the dependence of autowave velocity on parameters in isolated conditions is rarely studied, 13,16 and the VA relationship in the GS model in isolated conditions has never been reported, to the best of our knowledge. Here, we study this relation by the proposed analytical-numerical approach which, as we suggest, can be further adapted for more complex autowave models.

DESCRIPTION OF THE MODEL AND THE
ANALYTICAL APPROACH 2.1. Model Assumptions. We consider the GS model (Scheme 1) under the following assumptions: 1 Geometry and boundaries: The system is one-dimensional and isolated (there is no volume influx of components), and boundary influence is absent or negligible. 2 Reaction kinetics: All reaction kinetics are subjected to the law of mass action. Inhibition rate of the autocatalyst is nonzero, so the shape of its concentration profile is pulse (see Figure 2A, thick line). 3 Initial conditions: Initially, the precursor is uniformly distributed throughout the space, and the autocatalyst is absent everywhere except the "seed" quantity localized near one of the boundaries (for Figure 2A, near the left boundary). 4 Steady-state regime: The system has enough time to transit to the steady travelling wave regime, and transient process is not considered in this study. 5 Diffusion: Precursor diffusion may be neglected as qualitatively unimportant process provided the ratio of diffusion coefficients of precursor and autocatalyst is not large. Actually, the main effect of precursor diffusion is some decrease of its concentration compared to the initial finite value in front of the autowave front (in contrast to the autocatalyst diffusion which enables medium activation and wave propagation). The effect of this simplification will be addressed at the end of this article (see the Discussion and Figure 4A,B). This assumption is a fortiori valid for systems with an immobilized precursor. 2.2. Governing Equations. The standard form of RD equations corresponding to the Scheme 1 in isolated spatiallydistributed conditions and at neglected diffusion of v, is with boundary conditions and initial conditions Here, u and v are concentrations of the autocatalyst and its precursor, k and h are the activation and inhibition rate constants, D is the diffusion coefficient of the autocatalyst, v 0 is the initial uniform precursor concentration, and u 0 (x) is the initial autocatalyst concentration, which is supposed to be nonzero only in the small but finite spatial domain. If both u and v diffuse, the equation for v in eq 1 should be appended by the D v ·v xx term. The nondimensional form of eq 1 can be obtained defining t= , and h̃= h/kv 0 2 . Omitting tildes gives (2) Boundary conditions remain the same, and δ 0 (x) ≈ 1 in the finite spatial domain and zero elsewhere. Here, F(u,v) = v·u 2 − h·u. Note that F u ′(u = 0) = −h < 0, which gives threshold properties to the system. 8,21 Typical solution of these equations is shown in Figure 2A. Finally, to study solutions travelling with a steady velocity, c, that is, stationary in the moving frame x̂= x̃− c·t̃so that (u,v)(x,t) = (u,v)(x), eq 2 can be rewritten as In this eigenvalue problem for c, prime (′) denotes the derivative with respect to x, and below x̂will be referred to as x.
2.3. VA Relation from the Homoclinic Orbit Approximation. In this article, we suggest the following approach to estimate the VA relation for autowave systems that are qualitatively similar to the GS model. Let us denote = ′ = p u u ( ) u x d d and m = max(u) in the solution of eq 3. In the (v,u,p) phase space, this solution is a heteroclinic orbit connecting steady states (1,0,0) and (v 2 ,0,0), corresponding to x → +∞ and x → −∞, respectively. Projection of the orbit to the (u,p) plane is a homoclinic orbit starting from (0,0) and returning to (0,0), see Figure 2B, for example. The homoclinic . If we construct the formula (ansatz) approximating the shape of the orbit and containing c and m, then after integration, we will obtain an approximate VA relation.
(a) At x → ±∞, we have: u → 0 and vu 2 ≪ hu. Thus, u approximately obeys the homogeneous equation not contain- The characteristic equation and its roots are Solution of eq 4 is Thus, the orbit at the origin has tangents (in the (u,p) plane, see straight lines in Figure 2B) (b) At 0 ≤ u ≤ m, we suggest that the ansatz approximating half-orbits at p > 0 and p < 0 has the form so that these halves smoothly join at the point (m,0). To do this, we need functions g i (u) with properties The last condition guaranties smoothness of the junction. In this study, we use the specific form of g i (u) that fulfills all these conditions At n = 2, this equation defines an ellipse arc ( Figure 1A). We assume that n i > 1, and that n i may depend on h. With the above g(u) form, the equation for v in eq 3 can be integrated Here, = y u m , n ≡ n 2 , Β and Γ are beta and gamma functions, I 1 = ∞, and I 2 = 1. Equation 11 uniquely binds n and I n ( Figure  1B, solid line). In practice, it is convenient (but not necessary)

ACS Omega
Article to use a simpler approximate formula. All I n values in this study are well below 10 (see Figure 3B, gray line), and for I n ∈ [1,10], the following relation takes place ( Figure 1B, dashed line). The value a ≈ 0.15 was obtained by numerical minimization of the integral of the absolute difference between I n given by eqs 11 and 12 for n ∈ [n*,2], where n* ≈ 1.093 is the numerical solution of the equation I n (eq 11) = 10.
(c) Another relation between v 1 and m can be found by considering the equation for u in eq 3 , their best fit by eqs 7 and 9 with λ 1 ≈ 0,054 and n 1 ≈ 5 at p > 0 and λ 2 ≈ −0,72 and n 2 ≈ 1.27 at p < 0 (gray solid lines), effect of ±5% variation of n 2 (dashed half-orbits), and tangents given by eq 6 with the same λ 1,2 (dash-dotted straight lines). (C) Numerical solutions of eq 2 at three h-values (solid lines). Analytical half-orbits (p < 0) given by eqs 7 and 9 with λ 2 (c,h) from eq 5 and n 2 (c,m,h) from eq 18 at the same h-values (dashed lines).

ACS Omega
Article Substituting u = m and g(m) = 0 gives p′p = 0 − λ 2 m· g(m) 2−n , which is 0 if n < 2. Thus, eq 13 becomes the simple equality of the inhibition and production rates This relation can be rewritten as

Reduced Nagumo
Approximation. Another way to analyze eq 3 is to reduce them to a simpler model which is unique in that it has exact analytical solution. This model is known as the cold (isothermal) flame model, 1,13 purely cubic autocatalysis model, 2,10 or reduced Nagumo model. 5,6 In the region of steepest gradients, the rate of v → u reaction prevails over those of u inhibition and diffusion, thus v + u ≈ 1. More accurately, this first integral can be obtained by summation equations for u and v in the limit D v /D → 1 and h → 0. Setting v ≈ 1 − u gives Equation 19 describes the travelling front connecting two stable steady states, 0 and u 3 . This equation can be solved explicitly by guessing that the phase trajectory is u′ = p(u) = A· u·(u − u 3 ). 1,2,5,7,10 Substitution of this guess into eq 19 gives that = A 3. RESULTS 3.1. Overall View of Numerical Solution. The representative numerical solution of eq 2 (for h = 0.04) is shown in Figure 2A. The autowave travels from left to right with steady velocity and shape. The v(x) profile has a switchoff type (thin curve), and the u(x) profile has a pulse type (thick Parameters of the bottom half-orbit λ 2 (left axis, gray) and n 2 (right axis, black). λ 2 and n 2 (closed squares and closed circles) from the best fittings of the numerically obtained orbits at simultaneously varying λ 2 , n 2 , and m. n 2 (open circles) from varying of this parameter solely (while λ 2 , eq 5, and m were taken from numerical results). λ 2 and n 2 (gray and black lines) calculated by eqs 5 and 18 from the numerically obtained c and m, which are plotted in Figure 4. (B) I n (h) and n 2 (h) (gray and black lines) found from numerically obtained c and m ( Figure 4) using eq 18. For comparison, n 2 found by fitting of orbits is also shown (closed circles). Rectangle shows h and I n bounds used in Figure 4B to plot the analytical VA relation.

ACS Omega
Article curve). The same data in the u−p plane are shown in Figure 2B (black solid line), and solutions at several h-values are also shown in Figure 2C (solid lines). All solutions plot nearly homoclinic orbits starting from (0,0) to the bottom quarterplane, intersecting the 0u axis from below and tending to (0,0) from the top quarter-plane. Ingoing of orbits to (0,0) ceases at low h ( Figure 2C) because of finiteness of spatial interval of numerical integration, but this region of wave tail does not influence wave behavior, and increase of spatial interval should make all orbits definitely homoclinic (not shown).

Comparison of Numerical Solution and
Analytical Orbit Approximation. Surprisingly, an ansatz constructed in eqs 7 and 9 gives very good approximation of numerical orbits for every h when parameters λ i , n i , and m are varied just to get the best fit. Gray solid lines in Figure 2B show fitting example for h = 0.04, and dash-dotted lines show tangents (eq 6). Parameters λ 2 and n 2 of these fittings at different h values are plotted in Figure 3A (closed squares and closed circles). Obtained λ 2 and m are close to those calculated from numerical solutions (see gray squares vs gray lines for λ 2 ). That is why similar fittings were obtained when only n 2 was varied but λ 2 and m were fixed to values obtained from numerical solutions ( Figure 3A, open vs closed circles). Thus, the constructed ansatz represents the shape of orbit well and may be used to study general properties of the system.
3.3. Sensitivity of the Orbit Shape to the Parameter n 2 . Equations 7 and 9 mean that the orbit trajectory at x → ±∞ and u → m is mostly linked to velocity (through λ i ) and amplitude, respectively. To check that n i influence is localized only in the intermediate zone, two dashed half-orbits in Figure  2B show the effect of ±5% change of the n 2 value. Indeed, the effect of variation of n 2 is maximal between points (0,0) and (m,0), in the region of the steepest u(x) gradient. However, there is no any "hypersensitivity", as 10% change of n 2 leads to ∼11% change of |p| max .
3.4. Ansatz Shape Parameter n 2 can be Estimated from the Numerical Velocity and Amplitude. The next aim is to obtain the ansatz shape parameter related to autowave velocity and amplitude using model equations (rather than by orbit fitting). Numerical values of c and m in the entire range of h values are plotted in Figure 4A (solid lines). Expectedly, the steady autowave solution only exists at h below some threshold value, h*. We have found that h* ≈ 0.07. At h > h*, the wave does not survive irrespective of the triggering stimulus, δ 0 (x), and c = m = 0 are plotted. The same data in the c−m coordinates are plotted in Figure 4B (bold solid line), considering h as a parameter. Because of threshold, the c−m line does not reach the origin. Assuming that the derived analytical VA relation (eq 16) may be applied to these data, I n (h) and n 2 (h) were found using eq 18 ( Figure 3B, gray and black solid lines). The principal difference between these n 2 (h) (black line) and n 2 found above by simply fitting the u−u′ halforbits (circles) is that we derived eqs 16−18 taking into account model equations for both u and v, and that we used the computed autowave velocity in eq 18. Half-orbits with parameters λ 2 and n 2 calculated from numerically obtained c and m (eqs 5 and 18) are plotted in Figure 2C (dashed lines). They do not match numerical orbits precisely but at x → +∞ and at u → m they do.
Obtained I n (h) and n 2 (h) dependences are sloping quite mildly except for the region h → 0 ( Figure 3B). Because at h = 0, the u(x) profile in solutions of model equations has a shape of travelling front but not pulse, below we only consider nonzero h values.
3.5. Analytical VA Relation Weakly Depends on Particular (h,I n (h)) Values. Figure 3B shows that at almost all h (except h → 0), the values of I n (h) and n 2 (h) are bounded in the narrow ranges. For example, 2.6 ≤ I n ≤ 3.7 and 1.24 ≤ n 2 ≤ 1.34 at 0.01 ≤ h ≤ 0.07 (shaded rectangle). Substitution of these h and I n bounds into eq 17 gives the narrow corridor bounded by the dotted and dash-dotted gray lines in Figure 4B, which contains almost all numerical c−m results (bold solid line). Thus, any line passing through this corridor is an appropriate approximation of the c−m relation in the selected  Figure 3B.

ACS Omega
Article range of h values. For example, solid gray line in Figure 4B is plotted for intermediate values (h,I n (h)) = (0.04,2.9). The deviation is small even for h → 0 (m → 1 in Figure 4B).

DISCUSSION
VA relation is an important but underestimated characteristic of any autowave system. Usually, velocity is the only autowave property under consideration: first, nonzero velocity shows the very possibility of signal propagation, and its value gives the time required for a signal to cover the distance; second, the biological system response is usually regarded as a Yes/No (All/Nothing) type, so the precise value of autowave amplitude looks somewhat unimportant. However, autowave amplitude can be of major interest in various situations, such as the following: 1 The main outcome of a system is determined by autowave amplitude. For example, in blood clotting, the thrombin autowave leaves behind the polymerized fibrin mesh. 22,23 The fate of clot (to reside competing with blood flow or to embolize partly or wholly) depends on its physical, that is, rheological, properties such as storage and loss shear moduli. These macroscale properties depend on microscale characteristics of fibrin mesh (fibers' diameters, lengths, density, branching density, and cross-linking), which are determined mainly by the local transient thrombin concentration. 24−26 2 Autowave amplitude has great value in the economic, humanistic, and so forth, sense. For example, it designates the number (or fraction) of infected persons during spatial spreading of a disease. 6,19 3 The amplitude is easier to measure then the velocity, or the velocity is too low/too large to measure it correctly. In this case, the VA relation (like eq 17 for the GS model) can be used to estimate the velocity from the measurements of amplitude. In this study, we propose an approach to investigate VA relation in autowave systems in steady travelling regime. The key distinction of this approach is that it does not regard derivation of separate formulas for velocity and amplitude as functions of model parameters. Instead, it presumes the following steps: 1 Approximating the orbit shape in the phase space by a "petal" or "arc" nonlinear function (ansatz), which contains wave velocity and amplitude as well as a free parameter of orbit curvature. 2 Substituting the ansatz into model equations and integration of equations all along the wave front to the wave maximum; after some algebra, this gives an analytical VA relation. 3 Calculating dependence of the curvature parameter of the ansatz on model parameters from the derived analytical and numerical VA relations. 4 Checking whether the curvature parameter of ansatz has little effect on VA relation and may be regarded as fixed. In our study, the free parameter, n 2 , and its function, I n (eq 11), turned out to be bounded within the narrow ranges for all values of the model parameter h, allowing travelling wave solution (i.e., except for h → 0), see rectangle in Figure 3B. Any (h,I n (h)) pair from this range can be used in the analytical VA formula (eq 17) to well match the numerical results for all h ( Figure 4B). This agreement in VA relation did not require perfect matching between analytical and numerical orbits (see Figure 2C), suggesting that the orbit shape is the "finer" property of the system compared to the more "robust" VA relation property.
Moreover, the orbit shape may be approximated in numerous ways, and the choice of the function g(u) looks arbitrary. Indeed, g(u) in eq 9 may be multiplied by, for example, 1 − αu 2 or exp(−αu 2 ), still fulfilling the required conditions (eq 8) but introducing an additional free parameter, α. For the model system considered in this study, this complication looks unnecessary. Another advantage of the chosen g(u) form (eq 9) is the integrability of I n (eq 11), although this is not strictly necessary for the method, as eqs 16 and 17 contain I n but not n 2 . The question of the preferred choice of g(u) functions at different kinetic schemes will be addressed in a separate study.
Ansatz construction is a powerful but unformalized approach to obtain approximate and exact solutions of nonlinear differential equations. 7 The idea is to "guess" the structure of a solution basing on the knowledge of its general properties, and to test this guess by substituting it into the model equations. For example, the generalized Kolmogorov− Petrovskii−Piskunov−Fisher (KPPF) equation, that is eq 19 with F = u·(1 − u q ), can be solved by constructing an ansatz u(x) inspired by the main term of the approximate solution of the "classical" KPPF equation having q = 1. 7 Another example is suggestion of the polynomial form for the heteroclinic phase trajectory, p(u), if the reaction term in eq 19 is polynomial. 7 For cubic kinetics, this means that the phase trajectory can be the second order parabola, which is often used as a "guess", allowing an exact solution. 5 However, for autocatalytic systems with two variables (Scheme 1), this approach can only be used in case of approximately equal diffusion coefficients and zero inhibition rate, allowing to reduce a model to one equation (like eq 19). 1,2,10,13 Thus, all that has been said above applies to the models which are reducible to one equation and have front-like solutions. So far, there has not been constructed an ansatz for a model with two variables and nonzero inhibition rate of the autocatalyst.
Cubic autocatalysis with linear inhibition is a well-known model system arising in numerous chemical and biological tasks, from flame propagation to glycolysis. Even in the homogeneous case, this model was shown to have "exotic" (for chemistry) properties such as bistability and time-periodic solutions, provided there exists a continuous influx of the first component, that is, the conditions are open. 12,14,15 After the last articles, this model is sometimes called the GS model. The spatially distributed GS model was studied in several articles, 9,17,18,20 again in open conditions. In short, Petrov et al. 17 studied regimes of autowave propagation, reflection from boundaries, and from each other in the presence of influx of both components; Davidson et al. 18 studied fungal mycelia growth at the continuous influx of v (replenishment of the substrate); Hale et al. 20 derived the formula for the orbit at the large influx rate of v, equal to (v 0 − v)hD v /D (in notations of eq 1); Muratov and Osipov 9 used the singular perturbation technique to estimate autowave amplitude and velocity at the nonzero influx/removal rate, k 0 ; and their solutions diverge at k 0 → 0. Equations similar to the GS model emerged in the susceptible-infective model with the nonlinear infection rate, 19 with the only difference that the precursor influx is replaced by the reproduction rate of susceptibles. However, in isolated conditions, the GS model was considered only in two available studies. 13,16 ACS Omega Article Novozhilov and Posvyanski theoretically studied the cold flame propagation experimentally observed earlier by Voronkov and Semenov in the 0.03% CS 2 burning in the excess of oxygen, outside the ignition peninsula. 1,13 Considering the branched chain reaction scheme, they derived RD equations similar to the eq 3 (but with both components diffusible) and solved them numerically. Merkin and Needham 16 studied the same system analytically and numerically. In both studies, solutions existed only for h < h**, and the c(h) dependences were similar to that obtained in this study ( Figure 4A, gray solid line). To make the direct comparison, the v xx term was added into the equation for v in eq 2. Results of numerical integration of the "extended" system are also shown in Figure 4A (dashed lines). Solution exists for h < h** ≈ 0.045, in agreement with h** = 0.04611 13 and 0.0465. 16 At the zero inhibition rate, velocity c(h = 0) ≈ 0.7 also agrees with the published results 13,16 and with the reduced Nagumo approximation, (eq 20, dotted line in Figure 4A). Comparison of solution with and without diffusion of v ( Figure 4A,B, dashed vs solid curves) shows that diffusion of v, at least in the range 0 ≤ D v /D ≤ 1, affects properties of the system quantitatively but not qualitatively: it reduces autowave velocity and amplitude as well as h* but does not disrupt V(h), A(h), and VA relations. Novozhilov and Posvyanski obtained similar V(h) curves when varied the ratio D/D v in the range 0.5−2. 13 Qualitative unimportance of diffusion of v (at not very high ratio D v /D) can be well understood as it acts to smooth the v(x) profile: to decrease v at the wave front (where v ≈ 1) and to increase v in the wave tail (where v ≈ 0). Some decrease of v compared to v ≈ 1 at the wave front can indeed decrease the velocity, but the wave tail does not influence the system behavior. On the contrary, diffusion of u is a qualitatively important process as it brings u to the region where u ≈ 0, and the medium is activated in front of the front (see Figure 2A, x ≳ 85). Relative (qualitative) unimportance of v diffusion is the main reason we neglected it in this study. Another reason is that it makes equation for v ordinary and thus easier to integrate. Similar simplifications are used in other systems with low-mobile precursors, such as disease spreading (for susceptibles), 6,19 predator-prey (for pray), 6 resourceconsumer (for resource), 7,27 and so on.

NUMERICAL METHODS
To integrate eq 2 numerically, we used the Stormer−Encke method for space discretization (uniform mesh with N = 401 nodes, zero Neumann boundary conditions at x = 0 and x = 100) and the embedded Runge−Kutta−Fehlberg method of order 2(3) with automatic step size control for integration in time (mixed local error estimation with max norm, tol = 0.01, fac = 0.8, facmax = 5). 28 Increasing N and decreasing tol had negligible influence on c and m.

Funding
The publication has been prepared with support of the "RUDN University Program 5-100".

Notes
The author declares no competing financial interest. ■ ABBREVIATIONS