Automatic Proper Orthogonal Block Decomposition method for network dynamical systems with multiple timescales

Commun Nonlinear Sci Numer Simulat 131 (2024) 107844 Available online 12 January 2024 1007-5704/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Communications in Nonlinear Science and Numerical Simulation journal homepage: www.elsevier.com/locate/cnsns Research paper Automatic Proper Orthogonal Block Decomposition method for network dynamical systems with multiple timescales A. Banderaa,b,∗, S. Fernández-Garcíaa,b, M. Gómez-Mármola, A. Vidalc aEcuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, Calle Tarfia s/n, Seville, 41012, Spain bIMUS, Universidad de Sevilla, Calle Tarfia s/n, Seville, 41012, Spain cLaboratoire de Mathématiques et Modélisation d’Évry (LAMME), Univ Evry, CNRS, Université Paris-Saclay, IBGBI, 23 Bld de France, Evry, 91037, France ARTICLE INFO Keywords: Slow–fast dynamics Coupled oscillators Synchronization Network neuron model Reduced order models Mixed-mode oscillations ABSTRACT In this work, we introduce a novel reduced order model technique, based on the Proper Orthogonal Decomposition method, for dynamical systems with multiple timescales. The main ideas are to retain the structure of the original model, which is lost in the original POD procedure, while producing a competitive reduction in the number of equations and computational time, and to determine the best structure for the reduced system automatically, via a data-driven analysis of the original model data. For these novel techniques, we present some numerical tests for various behaviors of three different neural network models with multiple timescales, which support the use of these new methods. 1. Introduction Extensive research over the past three decades have been revolutionary by conceptualizing, analyzing and simulating network models of neuronal activity. Concerning the sole formalism of dynamical systems coupled as networks to represent interactions between assemblies of neurons, various paradigms and scales of modeling have been introduced depending on the questioning to be tackled, such as neural mass [1–3], neural fields [4,5], mean-fields [6] and dendritic trees [7], among others. However, common and still unsolved challenges arise in these approaches, especially the difficulty to classify and analyze the features of cellular activity synchronization, and the cost of reliable simulations for realistically sized networks. The first challenge mentioned above mainly implies results from the different timescales involved in neuronal activities, which are naturally translated into the models by the use of multiple timescale systems (for reproducing the individual cellular activities) coupled as a network. Synchronization mechanisms of multiple timescale systems are quite specific and different from synchronization of harmonic oscillators (especially the impact of the coupling types and strengths on the synchronization). Nevertheless, it is well-known that synchronization of neuron assemblies plays a key role in the neuronal information treatment and diffusion, and overall in the neural functioning. Several works (for instance, [8–17]) have addressed the question of analyzing synchronization mechanisms in the presence of canard phenomena and bursting oscillations, symptomatic of multiple timescale dynamics, in coupled oscillators. But the tremendously large panel of synchronization types and associated conditions remains to be analyzed, especially for large networks, as well as its physio-pathological interpretation. With this aim, efficient numerical techniques for producing reliable simulations are key tools and need to be thoroughly tested, due to the difficulty to catch structuring dynamical mechanisms specific to multiple timescale systems. As mentioned as a second ∗Corresponding author at: IMUS, Universidad de Sevilla, Calle Tarfia s/n, Seville, 41012, Spain. E-mail addresses: [email protected] (A. Bandera), [email protected] (S. Fernández-García), [email protected] (M. Gómez-Mármol), [email protected] (A. Vidal). https://doi.org/10.1016/j.cnsns.2024.107844 Received 1 September 2023; Received in revised form 22 December 2023; Accepted 8 January 2024 Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 2 A. Bandera et al. challenge above, simulation of network models might be computationally expensive when the number of cells is large. Therefore, intrusive Reduced Order Model (ROM) approximations have been developed [18–28]. For a state-of-the-art summary, we recommend Refs. [29,30]. These techniques aim to obtain an equivalent dynamical system of lower dimension that is a good approximation of the original Full Order Model (FOM), namely, the reduced problem. Despite the ability of these intrusive methods to approximate closely the solution generated by the original system [31], performing these techniques may make us lose the original structure of the problem, such as, for instance, the slow–fast separation of variables typical in neuroscience problems. This is due to the fact that the usual procedures to obtain the ROM do not consider any separation of variables, nor any given model structure. Therefore, we could be mixing fast and slow behaviors in the reduced model, which is not very natural. This problematic could also impede the direct study of the behaviors in the reduced framework. In order to preserve the structure of the problem, different strategies have been proposed in the literature. For instance, a technique based on Krylov subspaces was developed in the case of electronic circuits [19–21], and a technique to decompose the domain in computational simulations of fluids was introduced in [32]. In this paper, we present a different approach based, on the one hand, on the Proper Orthogonal Decomposition (POD) [18,22–27] and, on the other hand, on a new data-driven procedure based on a collinearity study to determine automatically the best structure for the reduced model, while obtaining a suitable reduction in computational time. The article outline is the following. In Section 2, we first summarize the POD method. We remark some properties, and we identify the issue of the structure loss. Then, we introduce in Section 2.1 a block approach to preserve the structure of the original problem after the reduction. We also adapt the properties of the POD method to this new block case. In Section 2.2, we present a data-driven procedure to determine the best block structure for the reduced model. After that, in Section 3we test and compare our methods for three neural network models with multiple timescales. In particular, Section 3.1 considers the case of a coupled network model of the Intracellular Calcium Concentration system, studied previously in [31]; Section 3.2 focuses on the study of a coupled network of Hindmarsh–Rose models [33–36] and Section 3.3 is devoted to the study of a coupled network model of Fold/Hopf bursters, found in pancreatic 𝛽-cells, among others [37–40]. We show the reduction of CPU time needed for the numerical integration of the reduced system compared to the original one, the sparsity indexes of the models as a measure of the preservation of the structure and an error comparison for the reduced models. Finally, Section 4is devoted to present conclusions and perspectives of the present work. Appendix is focused on analyzing the limits of these ROM techniques in some test cases. 2. Reduced order methods based on the proper orthogonal decomposition In this section, we begin by introducing the POD method. First, in order to build the ROM problem, we need to identify the number of equations. To do that, we can apply the POD method, one of the most popular approaches for reduced problem formulation [18,22,24–27]. Then, once we have the dimension of the reduced problem, the POD method also provides a basis for it, so we can formulate a ROM problem equivalent to the FOM problem. The POD method is based on the Singular Value Decomposition technique (SVD). It is known that the following property holds for the SVD of a matrix, S= (𝒘1|…|𝒘𝑛𝑠), see [41]: Proposition 2.1. Let 𝑁𝑟= {W∈R𝑁×𝑁𝑟∶W𝑇W=I𝑁𝑟}be the set of all orthonormal bases of 𝑁𝑟-dimensional subspaces of R𝑁. Then, 𝑛𝑠 ∑ 𝑖=1 ‖𝒘𝑖−UU𝑇𝒘𝑖‖2 2= min W∈𝑁𝑟 𝑛𝑠 ∑ 𝑖=1 ‖𝒘𝑖−WW𝑇𝒘𝑖‖2 2= 𝑟 ∑ 𝑖=𝑁𝑟+1 𝜎2 𝑖, where U∈R𝑁×𝑁𝑟arises from the SVD of Sand is a full matrix containing the first 𝑁𝑟left singular vectors in columns corresponding to the 𝑁𝑟greatest singular values and 𝑟is the number of non-zero singular values. Therefore, the SVD provides a basis that minimizes the sum of the square of the errors between each snapshot vector 𝒘𝑖and its projection onto the reduced subspace, and the error is equal to the sum of the squares of the neglected singular values. From now on, we focus on POD applied to dynamical systems. For a general parametric dynamical system  𝒙=𝑓(𝒙,𝜸), with 𝜸 the parameter vector, the POD method can be summarized as follows: 1. To solve the original problem for a certain partition of the parameter space, where the solution is qualitatively similar, and store the solutions in time in a matrix S, the so-called snapshot matrix. 2. To compute the SVD of the snapshot matrix S, and sort, in descending order, the singular values and its corresponding left singular vectors. 3. To apply Proposition 2.1 as a measure of the information lost to select the least amount of basis functions that retain the most information, according to a certain fixed tolerance. This will provide a submatrix of U, the matrix of left singular vectors, that contains the basis functions, which we also denote by U. 4. To construct the reduced problem as  𝒛=U𝑇𝑓(U𝒛, 𝑡), where 𝒛=U𝑇𝒙. To recover the solution in the original space, just a left-multiplication by Uof the reduced solution is needed. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 3 A. Bandera et al. However, as said in the introduction, performing this technique may make us to lose the original structure of the problem. This is due to the fact that in the original problem we work with quasi-sparse matrices and vectors, but in the reduction this quasi-sparsity is lost, and we have full matrices and vectors, as the corresponding basis matrix Uis full. Let us provide an illustrative example. In the following, in order to illustrate the main features of the methods proposed in a clear way, we will consider a dynamical system of 𝑁∈Nequations, that could be separated into two groups. However, this analysis can be generalized to 𝑀 < 𝑁 groups. Without loss of generality, we can consider that the variables can be sorted by group. A dynamical system in this framework could read, ( 𝒙𝟏  𝒙𝟐)=(A11 A12 A21 A22 )( 𝒙𝟏 𝒙𝟐)+(𝒃𝟏 𝒃𝟐),(1) where 𝒙𝒊and 𝒃𝒊,𝑖= 1,2, are vectors of dimensions 𝑁1and 𝑁2, respectively, such that 𝑁=𝑁1+𝑁2and A𝑖𝑗 ,𝑖, 𝑗 = 1,2, are real matrices with consistent dimensions. For the sake of clarity, we will consider A11 and A22 as the identity matrices, A21 = 0 and A12 a full matrix, so (A11 A12 A21 A22 )=(I A12 O I ). If we consider system (1) and we let Ube the matrix arising from the POD of some snapshots of the system, then the reduced order system reads,  𝒛= A𝒛+𝒃𝒓,(2) where  A=(I+U𝑇 11A11U21 U𝑇 11A12U22 U𝑇 12A12U21 I+U𝑇 12A12U22 ),with U=(U11 U12 U21 U22 ),and 𝒃𝒓=U𝑇(𝒃𝟏 𝒃𝟐), and we have lost the structural information of the original problem and the possible separation between variables. To solve this issue, we present the Proper Orthogonal Block Decomposition (POBD) method subsequently. 2.1. Proper orthogonal block decomposition In order to maintain the main properties of the FOM system when we perform the reduction, we have to take care of the structure of the matrix U. Taking into account the separation of variables in the FOM system (1), we aim to build a reduced basis for each set of variables. With that in mind, we divide the snapshot matrix S= (𝒘1|…|𝒘𝑛𝑠)into two blocks as follows S𝑇=(S𝑇 1|S𝑇 2). Now, we perform a SVD for each block, and we obtain a POD basis for each block, S1≈ 𝑟1 ∑ 𝑖=1 𝜎(1) 𝑖𝒖(1) 𝑖(𝒗(1) 𝑖)𝑇,S2≈ 𝑟2 ∑ 𝑖=1 𝜎(2) 𝑖𝒖(2) 𝑖(𝒗(2) 𝑖)𝑇, where 𝒖(⋅) 𝑖is the 𝑖th left singular vector of the corresponding matrix. Proposition 2.1 holds for each block, so as in the original POD method, we have a suitable criterion for the selection of the reduced basis dimension for each variable. We name the reduced dimensions 𝑁𝑟1, 𝑁𝑟2. Now, we focus on the reconstruction of the original snapshot matrix as a whole. We note that each block matrix has the correct number of columns, so no change in the block right singular vectors 𝒗(⋅) 𝑖is needed. However, each block left singular vector 𝒖(1) 𝑖lays in R𝑁1so we need to complete them in order to have a vector in R𝑁as we must, and respectively, for vectors 𝒖(2) 𝑖. A direct concatenation of the block left singular vectors does not work, as we do not know how the product 𝒖(𝑘) 𝑖(𝒗(𝑙) 𝑖)𝑇behaves, and we have no control over it. A successful approach is to add as many zeros as needed to each block left singular vector, taking into account the position of the block matrix in the original matrix, i.e., 𝒖(1) 𝑖↦(𝒖(1) 𝑖 𝟎𝑁2),𝒖(2) 𝑖↦(𝟎𝑁1 𝒖(2) 𝑖). With this decomposition, we can recover the original matrix as S≈(U1O𝑁1,𝑁𝑟2 O𝑁2,𝑁𝑟1U2)( 𝛴1O𝑁1,𝑁𝑟2 O𝑁2,𝑁𝑟1𝛴2)( V𝑇 1 V𝑇 2)=(S1 S2). If we select the matrix U=(U1O𝑁1,𝑁𝑟2 O𝑁2,𝑁𝑟1U2),(3) as the new basis for the reduced subspace, we are able to maintain the properties of the FOM system (1), as follows:  A=U𝑇AU =(U𝑇 1O𝑇 𝑁2,𝑁𝑟1 O𝑇 𝑁1,𝑁𝑟2 U𝑇 2)(I𝑁1A12 O𝑁2,𝑁1I𝑁2)( U1O𝑁1,𝑁𝑟2 O𝑁2,𝑁𝑟1U2) =(I𝑁𝑟1U𝑇 1A12U2 O𝑁𝑟2,𝑁𝑟1I𝑁𝑟2)=(I𝑁𝑟1  A12 O𝑁𝑟2,𝑁𝑟1I𝑁𝑟2). Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 4 A. Bandera et al. Table 1 Comparison of the sparsity index for between the FOM, POD and POBD systems. FOM Non-zero elements ≈𝑁1+𝑁2+𝑁1𝑁2 Total elements (𝑁1+𝑁2)2 Sparsity index ≈𝑁1(𝑁1− 1) + 𝑁2(𝑁2− 1) + 𝑁1𝑁2 (𝑁1+𝑁2)2 POD POBD Non-zero elements ≈𝑁2 𝑟≈𝑁𝑟1+𝑁𝑟2+𝑁𝑟1𝑁𝑟2 Total elements 𝑁2 𝑟(𝑁𝑟1+𝑁𝑟2)2 Sparsity index ≈ 0 ≈ 𝑁𝑟1(𝑁𝑟1− 1) + 𝑁𝑟2(𝑁𝑟2− 1) + 𝑁𝑟1𝑁𝑟2 (𝑁𝑟1+𝑁𝑟2)2 We remark that we keep the original structure for our reduced system. We name this method the Proper Orthogonal Block Decomposition (POBD). In order to measure the good behavior of this method, taking into account just the number of equations of the reduced system could be not enough, as we could also need to retain some information about the structure of the original matrix. Bearing that in mind, we propose a new index, that we name the sparsity index, defined as the quotient between the number of zero entries and the number of total elements of the matrix of a system of differential equations. Table 1 shows the comparison between the sparsity index of the FOM, POD and POBD systems for the academic test case (1). We can easily check that, for this test case, the sparsity index of the POBD systems is much more in concordance with the one of the FOM system, while the sparsity index of the POD system is close to zero. For the original POD method, we are able to obtain an error bound depending on the values of the neglected eigenvalues, see Proposition 2.1. In our case, we are able to obtain a similar error bound applying similar techniques. Proposition 2.2. Let Ube an orthonormal matrix built as in (3) for a dynamical system with 𝑁𝑏≥2blocks. Then, the next error bound holds 𝑛𝑠 ∑ 𝑖=1 ‖𝒘𝑖−U𝑇U𝒘𝑖‖2 2= 𝑁𝑏 ∑ 𝑗=1 𝑟𝑗 ∑ 𝑖=𝑁𝑟𝑗+1 (𝜎(𝑗) 𝑖)2. Proof. This proof is an adaptation of the proof of Proposition 2.1 in [41] to our case. We detail it for 𝑁𝑏= 2. For a matrix, A∈R𝑁×𝑛𝑠we have ‖A‖2 𝐹=∑𝑛𝑠 𝑖=1 ‖𝒂𝑖‖2 2, where 𝒂𝑖denotes the 𝑖th column of A. Thus, in our case we can write 𝑛𝑠 ∑ 𝑖=1 ‖𝒘𝑖−UU𝑇𝒘𝑖‖2 2=‖S−UU𝑇S‖2 𝐹.(4) Now, as S=(S1 S2)=(S1 O𝑁2,𝑛𝑠)+(O𝑁1,𝑛𝑠 S2), and UU𝑇S=(U1U𝑇 1S1 U2U𝑇 2S2)=(U1U𝑇 1S1 O𝑁2,𝑛𝑠)+(O𝑁1,𝑛𝑠 U2U𝑇 2S2), along with the facts that ‖𝑆‖𝐹=‖𝑆𝑇‖𝐹and ‖(𝐴|𝐵)‖2 𝐹=‖𝐴‖2 𝐹+‖𝐵‖2 𝐹, we can decompose the last term of (4) as ‖S−UU𝑇S‖2 𝐹=‖S1−U1U𝑇 1S1‖2 𝐹+‖S2−U2U𝑇 2S2‖2 𝐹.(5) Finally, applying Proposition 2.1 to each term, we obtain the desired result. We can remark that this procedure is suitable for any number of blocks and that when considering just 1 block, we recover the classic POD method. Remark 2.1. We remark that despite the similarities between the error bounds of Propositions 2.1 and 2.2, for a number of blocks 𝑁𝑏≥2the latter might not be the optimal for reducing onto a subspace of dimension ∑𝑁𝑏 𝑗=1 𝑁𝑟𝑗. 2.2. Automatic proper orthogonal block decomposition In problems with 𝑁𝑏large, there could be blocks of variables that could be studied as a whole block, allowing us to improve the error estimation presented in Proposition 2.2, as we would need to deal with a lesser number of blocks. In this section, we propose Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 5 A. Bandera et al. a procedure to identify those blocks that could be merged into a larger block by means of a collinearity analysis of the snapshots. We would name this method Automatic Proper Orthogonal Block Decomposition (APOBD). The procedure of the selection of blocks is summarized in the following. We have illustrated the steps with the academic example of the previous Eq. (1), 1. The first step is to compute the correlation matrix of the snapshots C∈R𝑁×𝑁. This could also be seen as the first step for the POD, as the singular values of a matrix are the square root of the eigenvalues of the respective correlation matrix, C=SS𝑇=(S1 S2)(S𝑇 1S𝑇 2)=(S1S𝑇 1S1S𝑇 2 S2S𝑇 1S2S𝑇 2)=(C11 C12 C𝑇 12 C22 ).(6) 2. Then, we compute the mean correlation of each block with itself and the others, so we end up with a symmetric matrix of dimension equal to the number of blocks. In order to check if the a priori blocks were well located, we can check that the elements of the diagonal of this matrix are close to 1. L=⎛⎜⎜⎜⎜⎜⎝ ∑𝑁1 𝑖,𝑗=1(C11)𝑖𝑗 𝑁2 1∑𝑁1 𝑖=1 ∑𝑁2 𝑗=1(C12)𝑖𝑗 𝑁1𝑁2 ∑𝑁1 𝑖=1 ∑𝑁2 𝑗=1(C12)𝑖𝑗 𝑁1𝑁2∑𝑁2 𝑖,𝑗=1(C22)𝑖𝑗 𝑁2 2 ⎞⎟⎟⎟⎟⎟⎠ =(𝑙11 𝑙12 𝑙12 𝑙22 ),with 𝑙𝑖𝑖 ≈ 1.(7) 3. As the third step, for each pair of different blocks (𝑖, 𝑗), we compute the collinearity 𝑐𝑙𝑖𝑗 , that is, the determinant of the 2 ×2 submatrix of the matrix Lof the previous step (7) corresponding to the blocks studied, 𝑐𝑙𝑖𝑗 = det (𝑙𝑖𝑖 𝑙𝑖𝑗 𝑙𝑖𝑗 𝑙𝑗𝑗 ).(8) If this value is lower than a certain fixed tolerance, we can consider the two blocks as a unique block. At most, we would need to compute 𝑁𝑏(𝑁𝑏− 1)∕2 coefficients as 𝑐𝑙𝑖𝑖 = 0 and 𝑐𝑙𝑖𝑗 =𝑐𝑙𝑗𝑖. 3. Test cases This section is devoted to the validation and numerical study of all the aforementioned methods, POD, POBD and APOBD. We will apply these reduction techniques to problems arising from the mathematical modeling of biological phenomena exhibited in cells or neurons: Intracellular Calcium Concentration model [31,42], Hindmarsh–Rose model [33–36] and pancreatic 𝛽-cells models [37,38,40,43]. We will work with network models with a high number of coupled oscillators, and we will compare the reduction techniques by means of some classical indicators, such as speed-up and number of reduced equations, and also through the sparsity index, an indicator of the reproduction of the structure of the original problems based on the sparsity of the reduced matrices associated to each reduced model. 3.1. Intracellular calcium concentration network model In this section, we test the aforementioned methods on the Intracellular Calcium Concentration model studied in [31,42]. In realistic problems, we have more than the two cells originally considered in those works. We consider here the case where we have 𝑁∈Ndifferent cells and symmetric coupling. We can model the system with 𝑁coupled oscillators as 𝑂𝑖 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑥𝑖=𝜏(−𝑦𝑖+𝑓(𝑥𝑖) − 𝜙𝑓(𝑧𝑖)), 𝑦𝑖=𝜏𝜀𝑘𝑖(𝑎0𝑥𝑖+𝑎1𝑦𝑖+𝑎2+1 𝑁∕2 𝑁 ∑ 𝑗=1 𝑐𝑖𝑗 (𝑥𝑖−𝑥𝑗)), 𝑧𝑖=𝜏𝜀 (𝜙𝑟(𝑥𝑖) − 𝑧𝑖−𝑧𝑏 𝜏𝑧), (9) where 𝑖= 1,2,…, 𝑁,𝑂𝑖denotes the 𝑖th oscillator formed by (𝑥𝑖, 𝑦𝑖, 𝑧𝑖)and 𝑐𝑖𝑗 ∈ [−1,1] with 𝑐𝑖𝑗 =𝑐𝑗𝑖 and 𝜙𝑓(𝑧) = 𝜇𝑧 𝑧+𝑧0 , 𝜙𝑟(𝑥) = 𝜆 1 + exp(−𝜌(𝑥−𝑥𝑜𝑛)) .(10) In the following section, for parameters that are fixed among the simulations, we have kept the values used in [31,42]. Therefore, we consider 𝑓(𝑥)=−𝑥3+ 4𝑥and 𝑎1= −0.1, 𝑎2= 0.8, 𝜀 = 0.06, 𝜏 = 37, 𝜇 = 2.4, 𝑧0= 5, 𝑧𝑏= 1, 𝜏𝑧= 2, 𝑥𝑜𝑛 = −0.45, 𝜆 = 1.75, 𝜌 = 4.5.(11) Thus, 𝑥would be the fast variable and 𝑦, 𝑧 the slow ones. In a first approach, we consider the cells separated following a modular network structure with two clusters, namely, first cluster, are the cells with indexes 𝑖= 1,…, 𝑁−, and second cluster, the cells with indexes 𝑖=𝑁−+ 1,…, 𝑁. This choice is based on the motoneurons network properties studied in the zebrafish embryo in [44], and has been considered previously in [31]. We set 𝑁−=⌊𝑁∕2⌋, that is, the greatest integer lower than or equal to 𝑁∕2, and 𝑁+=⌈𝑁∕2⌉, that is, the lowest integer greater than or equal to 𝑁∕2. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 6 A. Bandera et al. 3.1.1. A priori study In this section, we build the algebraic form of the original network problem (9). We denote 𝒘= (𝑥1,…, 𝑥𝑁, 𝑦1,…, 𝑦𝑁, 𝑧1,…, 𝑧𝑁), and we can write problem (9) as  𝒘=𝑓(𝒘,𝜸),with 𝒘(0) = 𝒘0,(12) where 𝜸= (𝑐𝛼, 𝑐𝛽,𝒌), with 𝑘𝑖the heterogeneity parameters, and 𝑐𝛼, 𝑐𝛽the intra-cluster and inter-cluster coupling, respectively. We note that we have changed the order of the variables in order to help the reader to understand the relationship between them, and that it has no consequences in the following study. First, for the sake of clarity, we introduce some notations. Notation 3.1. We denote by 𝟎𝑁the null column vector with 𝑁components, 𝟏the column vector with 𝑁components, and all its components 1 and 𝒌 the column vector whose entries are the heterogeneity parameters. On the other hand, I𝑁is the identity matrix of order 𝑁,O𝑁is the null matrix of order 𝑁and K= diag(𝒌)is a diagonal matrix whose entries are the heterogeneity parameters. Also, we define the following parameters, 𝛴𝑖={(𝑁−− 1)𝑐𝛼+𝑁+𝑐𝛽,for 𝑖= 1,…, 𝑁−, (𝑁+− 1)𝑐𝛼+𝑁−𝑐𝛽,for 𝑖=𝑁−+ 1,…, 𝑁. Note that, when 𝑁is even, 𝛴𝑖=𝛴= (𝑁∕2 − 1)𝑐𝛼+𝑁𝑐𝛽∕2,∀𝑖= 1,…, 𝑁. Finally, we consider, M(𝑐𝛼, 𝑐𝛽) = ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 𝛴1−𝑐𝛼… −𝑐𝛼−𝑐𝛽… … −𝑐𝛽 −𝑐𝛼⋱ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋱ −𝑐𝛼⋮ ⋮ −𝑐𝛼… −𝑐𝛼𝛴𝑁−−𝑐𝛽… … −𝑐𝛽 −𝑐𝛽… … −𝑐𝛽𝛴𝑁−+1 −𝑐𝛼… −𝑐𝛼 ⋮ ⋮ −𝑐𝛼⋱ ⋱ ⋮ ⋮ ⋮ ⋮ ⋱ ⋱ −𝑐𝛼 −𝑐𝛽… … −𝑐𝛽−𝑐𝛼… −𝑐𝛼𝛴𝑁 ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ Thus, top left submatrix is 𝑁−×𝑁−, top right is 𝑁−×𝑁+, bottom left is 𝑁+×𝑁−, and bottom right is 𝑁+×𝑁+. We remark here that the matrix M(𝑐𝛼, 𝑐𝛽)admits an affine decomposition depending on 𝑐𝛼, 𝑐𝛽and 𝛴𝑖. Injecting the previous notation in Eq. (12), we can rewrite system (9) as  𝒘=𝜏((A𝑛𝑝 +A𝑝(𝒌) + B(𝑐𝛼, 𝑐𝛽,𝒌))𝒘+𝒃𝑛𝑝 +𝒃𝑝(𝒌) + 𝒃𝑛𝑙(𝒘)),(13) where A𝑛𝑝 and 𝒃𝑛𝑝 represent the non-parametric part, A𝑝(𝒌)and 𝒃𝑝(𝒌)the parametric part depending on the heterogeneity parameters, B(𝑐𝛼, 𝑐𝛽)is the parametric part relative to the coupling and 𝒃𝑛𝑙(𝒘)denotes the nonlinear part of system (9). They can be expressed as follows A𝑛𝑝 =⎛⎜⎜⎝ 4I𝑁−I𝑁O𝑁 O𝑁O𝑁O𝑁 O𝑁O𝑁−𝜀∕𝜏𝑧I𝑁⎞⎟⎟⎠ ,A𝑝(𝒌) = 𝜀⎛⎜⎜⎝ O𝑁O𝑁O𝑁 𝑎0K𝑎1K O𝑁 O𝑁O𝑁O𝑁⎞⎟⎟⎠ , B(𝑐𝛼, 𝑐𝛽) = 𝜀 𝑁∕2 ⎛⎜⎜⎝ O𝑁O𝑁O𝑁 KM(𝑐𝛼, 𝑐𝛽)O𝑁O𝑁 O𝑁O𝑁O𝑁⎞⎟⎟⎠ , 𝒃𝑛𝑝 =𝜀𝑧𝑏 𝜏𝑧⎛⎜⎜⎝ 𝟎𝑁 𝟎𝑁 𝟏𝑁⎞⎟⎟⎠ ,𝒃𝑝(𝒌) = 𝜀⎛⎜⎜⎝ 𝟎𝑁 𝑎2𝒌 𝟎𝑁⎞⎟⎟⎠ ,𝒃𝑛𝑙(𝒘) = ⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝ −𝑥3 1−𝜙𝑓(𝑧1) ⋮ −𝑥3 𝑁−𝜙𝑓(𝑧𝑁) 𝟎𝑁 𝜀𝜙𝑟(𝑥1) ⋮ 𝜀𝜙𝑟(𝑥𝑁) ⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠ . With a first look into the model, we can easily see that system (13) depends on many parameters. In the following study, we set random heterogeneity parameters 𝑘𝑖∈ [1,1.5] following a truncated normal distribution (1.25,0.25), fixed for all the simulations, and fix 𝑐𝛼= 1, that is, we consider that the cells within each cluster are coupled in an excitatory fashion. Therefore, we perform our reduction with respect to the inter-cluster coupling parameter 𝑐𝛽, that is, the coupling between cells of different clusters. This selection is made to be in concordance with previous works [31], and to avoid chaos and bizarre behaviors. As a matter of fact, chaotic behaviors are tremendously sensitive to perturbations. Consequently, a test case for any reduction process aiming at accelerating the simulations cannot be successful in the presence of chaos. Nevertheless, we comment this case and show corresponding simulations in Appendix A.1. Before proceeding to the computational study, we need to set which behaviors we are trying to model, along with its range on the parameter 𝑐𝛽and also the number of blocks for the POBD. We will focus on three main behaviors: Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 7 A. Bandera et al. Fig. 1. Representation of the different behaviors to study. In blue, the cells belonging to the first cluster and in red, the cells belonging to the second cluster. We only represent 𝑥and 𝑧, for the sake of clarity. 1. Almost-in-Phase synchronization (AiP). For 𝑐𝛽∈ [0.1,1]. In this range of the parameter, each cell reaches asymptotically an orbit similar to the cycle existing for the relaxation oscillation of the uncoupled case without Mixed Mode Oscillations (MMOs). All the 3D oscillators are synchronized with a small phase-shift or change of frequency depending on its heterogeneity parameter. This behavior is shown in Fig. 1(a). 2. Antiphase Synchronization (AS). For 𝑐𝛽∈ [−0.45,−0.05]. In this range of the parameter, each cell within each cluster reach asymptotically the same orbit, and the phase-shift between the clusters is close to half the period with a small phase-shift or change of frequency depending on its heterogeneity parameter. This behavior is shown in Fig. 1(b). 3. Total Oscillation Death (TOD). For 𝑐𝛽∈ [−1,−0.55]. In this range of the parameter, each cell within each cluster reaches asymptotically a stable equilibrium. This behavior is shown in Fig. 1(c). As we have two clusters with three types of variables we will perform the POBD with 3 and 6 blocks, the first POBD will have one block per variable (𝒙,𝒚,𝒛), and we will name it POBD3, and the second POBD will have one block per cluster variable (𝒙𝐼,𝒙𝐼𝐼 ,𝒚𝐼,𝒚𝐼𝐼 ,𝒛𝐼,𝒛𝐼𝐼 ), where subindex 𝐼denotes the variables relatives to the first cluster and subindex 𝐼𝐼 denotes the variables relatives to the second cluster. We will name this method POBD6. As APOBD groups the blocks automatically, we start with the maximum number of possible blocks, in our case 6, one per cluster and variable. 3.1.2. Numerical results In this section, we will perform various numerical tests for different number of cells, namely, 𝑁= 200,500,1000,2000. In every test case, we will let the original system (13) compute from 𝑇= 0 to 𝑇= 50 via a fourth order Runge–Kutta solver with fixed time step of ℎ= 10−3. We perform every method with the data from 𝑇= 35 to 𝑇= 50 as we need that the system is almost at an asymptotic behavior. In the construction of reduced systems, the tolerance criteria for the collinearity analysis in the APOBD method is set at 0.2, that give us two blocks in AiP behavior, one for the fast variables (𝒙)and another one for the slow ones (𝒚,𝒛), four blocks in AS behavior, (𝒙𝐼),(𝒙𝐼𝐼 ),(𝒚𝐼,𝒛𝐼)and (𝒚𝐼𝐼 ,𝒛𝐼𝐼 ), that is, one block per timescale and cluster and just one block for TOD behavior, that is, all the variables mixed together. The tolerance criteria of the selection of the reduced dimension in POD methods is set to 10−9for AiP behavior and 10−6for AS and TOD behaviors. Nonlinear terms of (13) are addressed through the Empirical Interpolation Method [45,46] (from now on, EIM) with a tolerance of 10−6. This is a technique that allows us to compute the nonlinearities of the system as a linear combination of some basis. In the following, we compare the number of equations, the sparsity of the reduced model via the sparsity index, and the computational time spent to solve the reduced system from 𝑇= 50 to 𝑇= 70, via the speed-ups, defined as the quotient between the computational time of the original system (13) and the computational time of the reduced system. The objectives of this test case are the following: First, to validate the new proposed methods. In order to do so, we compare all the methods with 𝑁= 200 cells in the three behaviors mentioned above. Then, in order to compare in more detail the methods, we will test them with a more large amount of cells, that is, 𝑁= 500,1000,2000 for AiP and AS behaviors, as the TOD behavior has not a major interest, as it reproduces the system in a steady state and the reduction is straightforward. In Table 2, we show the simulation data of the original system (13) for every test case. Validation of the methods for 𝑁= 200 The results of the reductions are summarized in Table 3. We see that for this test, POD performs better in terms of speed-up just for AiP behavior, but in terms of sparsity POBD6 provides the best index in the three cases, that is, the best separation of variables. We see that both POBD3 and APOBD provide a well trade-off between sparsity and speed-up. For AS behavior, POBD6 and APOBD Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 8 A. Bandera et al. Table 2 Data for system (13) for different test cases. Test case 𝑁= 200 𝑁= 500 𝑁= 1000 𝑁= 2000 Number of equations 600 1500 3000 6000 Sparsity index 0.8867 0.8880 0.8884 0.8887 CPU time ≈6.2 s ≈14 s ≈126 s ≈475 s Table 3 Comparison of the data for the reduced system for the original system (13) with 𝑁= 200. Notation a+...+c indicates the number of equations in each block. In bold fashion, the best results for each behavior. Behavior Method POD POBD3 POBD6 APOBD Number of equations 33 26+12+10 21+23+10+11+8+9 26+16 AiP Sparsity index 0 0.6510 0.7518 0.3685 Speed-up 6.04 4.20 2.70 4.24 Error (⋅10−3)5.06 4.46 6.95 4.35 Number of equations 57 53+14+8 27+26+7+7+4+4 27+26+9+8 AS Sparsity index 0 0.6874 0.7739 0.6937 Speed-up 2.02 1.67 2.20 2.17 Error (⋅10−3)3.59 7.80 6.24 6.98 Number of equations 1 1+1+1 1+1+1+1+1+1 1 TOD Sparsity index 0 0.4444 0.6667 0 Speed-up 34 15.32 9.29 40.06 Error (⋅10−4)3.09 0.40 0.21 3.09 Table 4 Comparison of the number of equations of the reduced systems for the original system (13) in AiP behavior. Number of cells POD POBD3 POBD6 APOBD 500 36 29+14+11 25+27+12+13+10+11 29+18 1000 39 32+14+12 28+31+12+14+10+12 32+18 2000 42 36+15+12 35+34+15+14+12+11 36+18 are the methods that provide the best speed-ups and sparsity indexes. For TOD behavior, APOBD provides the best speed-up, while POBD6 provides the best separation of variables and error. For illustration purposes, in Fig. 2 we show the first 3 basis obtained for each method and each block in AS behavior. Note that the basis for POD (Fig. 2(a)) and POBD3 (Fig. 2(b)) are much more erratic, and we cannot identify neither the separation between variables nor clusters, while for POBD6 (Fig. 2(c)) and APOBD (Fig. 2(d)) this separation is clearer, due to the higher sparsity index. This means that, we could not need to go back to the full-dimensional space to analyze the data, allowing us to save more computational time. At this point, we can affirm that this validation test is finished successfully. The reduced methods proposed are able to reproduce the solution of the original problem while performing a reduction in computational times. It should be noted that while all the methods are comparable in terms of sparsity and error, the novel methods also provide us with a new feature, the perseverance of the structure of the original system. Now, we are able to continue our study. In next sections, we will focus just on the AiP and AS behaviors. As mentioned in the previous, the TOD behavior is not of significant interest as it simply replicates the system in a steady state, and the reduction process is straightforward. Comparison for larger number of cells in AiP behavior In this section, we study in detail the reduction in the AiP behavior of system (13) for larger numbers of cells, 𝑁= 500,1000,2000. The results are summarized in Table 4, which shows the number of variables of each cell number case, and Figs. 3(a)–3(c), which represent the evolution of sparsity, speed-up and error, respectively. From Table 4, we can see that APOBD minimizes the number of equations with respect to the POBD methods, while maintaining the structure of the system. The best method in terms of sparsity is POBD6 in every case, see Fig. 3(a), however, a better sparsity index does not imply a better speed-up, as it is also the worst method in terms of speed-up, see Fig. 3(b). It seems that the best method in terms of the tradeoff between sparsity and speed-up is APOBD. We can observe that the relative errors are comparable, as they have the same order, see Fig. 3(c). Comparison for larger number of cells in AS behavior In this section, we study in detail the reduction in the AS behavior of system (13) for larger numbers of cells, 𝑁= 500,1000,2000. The results are summarized in Table 5, which shows the number of variables of each cell number case, and Figs. 4(a)–4(c), which represent the evolution of sparsity, speed-up and error, respectively. From Table 5, we can see again that APOBD minimizes the number of equations with respect to the POBD methods, while maintaining the structure of the system. The best method in terms of sparsity is POBD6 in every case, see Fig. 4(a), however, a better sparsity index does not imply a better speed-up, as it is not the best method in terms of speed-up, see Fig. 4(b). We can see that methods POBD3 and APOBD Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 9 A. Bandera et al. Fig. 2. Comparison of the basis for each method. In all the cases, the first basis is in blue continuous fashion, the second in red dash-dotted fashion and the third basis in yellow dashed fashion. (a) POD. (b) POBD3. (c) POBD6. (d) APOBD. Fig. 3. Comparisons of the results for each method in AiP behavior of system (13). provide a similar sparsity index, but APOBD is much better in terms of speed-up. Again, it seems that the best method in terms of the tradeoff between sparsity and speed-up is APOBD. Again, the relative errors are comparable between the methods, see Fig. 4(c). It is clear that our proposed method APOBD gives the best results in terms of the trade-off between model reduction, saving of computational time and explainability of the reduced variables, as it needs a number of equations that is very low in comparison Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 16 A. Bandera et al. Fig. 11. FHB behavior of system (18) for 𝑁= 100 cells. Table 9 Simulation data for system (18). Test case 𝑁= 100 𝑁= 500 𝑁= 1000 Number of equations 300 1500 3000 Sparsity index 0.8833 0.8867 0.8887 CPU time ≈4 s ≈65 s ≈280 s Table 10 Comparison of the number of equations of the reduced systems for the original system (18) in FHB behavior. Number of cells POD APOBD 100 3 2+2+2 500 3 2+2+2 1000 3 2+2+2 3.3.2. Numerical results In the following study, we will consider various instances for the number of cells 𝑁= 100,500,1000. In every test case, we will let the original system (18) compute from 𝑇= 0 to 𝑇= 400 via a fourth order Runge–Kutta solver with fixed time step of ℎ= 10−2. And we perform both reduction methods with the data from 𝑇= 200 to 𝑇= 400 as we need that the system is almost at an asymptotic behavior, this will be enough to capture the main behavior of the system as shown in Fig. 11. In the construction of reduced systems, the tolerance criteria for the collinearity analysis in the APOBD method is set at 0.2, that gives us three blocks in the behavior of the study. The tolerance criteria of the selection of the reduced dimension in POD methods is set to 10−6for every behavior, and nonlinear parts of (18) are addressed by means of the EIM with a tolerance of 10−6. Table 9 summarizes the simulation data of the original system (18) for the test case. Now, we compare the number of equations, the sparsity index, and the speed-ups, from 𝑇= 400 to 𝑇= 800, for the original model and the reduced ones. The results are summarized in Table 10, which shows the number of variables of each cell number case, and Figs. 12(a)–12(c), which represent the evolution of sparsity, speed-up and error, respectively. From Table 10, we can see again that the reduction is stronger when the number of original equations rises, as the number of reduced equations remains the same. We can again see, that although the speed-up is better for the POD method (see Fig. 12(b)), APOBD provides three additional orders in the relative error, Fig. 12(c). So, for almost the same computational time, we obtain a much more precise approximation of the solution of the original problem (18). For illustration purposes, in Fig. 13 we show the first basis obtained for each method and each block in the studied behavior, in the three basis for POD (Fig. 13(a)), where we cannot see a separation in the timescales, and in the basis for APOBD (Fig. 13(b)), where this separation is clearer, thanks to the higher sparsity index. 4. Conclusions and perspectives In this work, we have developed new ROM techniques based on the POD method for dynamical systems with multiple timescales. They are able to maintain the structure of the original model, while producing a competitive speed-up. In this case, we have Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 17 A. Bandera et al. Fig. 12. Comparisons of the results for each method in FHB behavior of system (18). Fig. 13. Comparisons of the basis for each method obtained in the studied behavior. In all the cases, the first basis is in blue continuous fashion, the second in red dash-dotted fashion and the third basis in yellow dashed fashion, this one only for the POD method. (a) POD. (b) APOBD. performed an a priori analysis of how the variables of the original model could be grouped and then apply a data-driven analysis based on the correlation between these groups of variables in order to reduce the number of blocks to consider, while maintaining the most part of the structure of the original problem. We have also provided an error estimator of the information missed when performing these methods, depending on the number of considered blocks. Furthermore, we have tested these methods for several biological neural models with multiple timescales that exhibit a rich variety of different qualitative behaviors. These new methods have been proven to produce an accurate approximation of the original solution in much less computational time, and we obtain results in speed-up and error that are better, or at least, as good as, the original POD procedure to build the reduced models. These new methods also provide a new feature, they keep the structure of the original model, which could be useful in the sense that we can see the different timescales in the new basis of the reduced model and this allows us to exploit the separation of the timescales. Furthermore, we have studied the evolution of the performance of these methods augmenting the number of cells in the original models, and they seem to perform better the more cells are considered. We have also noticed that these new techniques seem to perform better when there are more fast variables than slow variables and where there is a switch between the fast and slow behaviors. Therefore, these new techniques based on the POD and a separation of variables to keep the structure of the original model are a useful tool for numerical study of network models, as they are able to reproduce complex behaviors that these models may exhibit. In future works, we aim to apply these new ROM techniques considering various biological parameters; we also intend to apply them to perform a numerical analysis of the repartition of the behaviors of the models presented, as well as to apply these new techniques to different problems arising from other fields, like physical problems and computational fluid dynamics problems. Furthermore, we intend to study how these reductions can be performed in the transitions between different behaviors. CRediT authorship contribution statement A. Bandera: Writing – original draft, Writing – review & editing. S. Fernández-García: Writing – original draft, Writing – review & editing. M. Gómez-Mármol: Writing – original draft, Writing – review & editing. A. Vidal: Writing – original draft, Writing – review & editing. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 18 A. Bandera et al. Fig. 14. Matrix Larising after computing the mean correlation of each block of variables of system (13) for 𝑁= 100 with 𝑐𝛼= −0.5and 𝑐𝛽∈ [−0.25,0.25]. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability No data was used for the research described in the article. Acknowledgment This work has been supported by the Spanish Government Project PID2021-123153OB-C21. Appendix. Limitations of the reduced order modeling techniques The methods presented in the previous sections have been proven to provide competitive results in terms of maintaining the structure of the original problem, studied via the sparsity index, providing a reduction in the computational time, studied via the speed-up and being true to the original solution, studied via the relative error. However, there are some limitations to these techniques, both older and novel. Here, we discuss the limitations of the methods presented in this article. A.1. Limitations in the network model of intracellular calcium concentration As said in the presentation of the APOBD method, we could check if the block were well located, just checking if the elements of the diagonal of the matrix were close to one. A natural question arises: What happens in the case where there is no correlation inside each block of variables? In order to illustrate the results in this case, we consider system (13) for 𝑁= 100 with 𝑐𝛼= −0.5, that is, inhibitory coupling within each cluster, and 𝑐𝛽∈ [−0.25,0.25]. For the study, we consider the same parameters test as the ones used for behavior AS, and we compare POD, POBD3, POBD6, APOBD. The APOBD pre-analysis provides the matrix Lshown in Fig. 14. It is easy to see that the elements of the diagonal are far from zero, so the POBD methods and the APOBD method, would not provide a faithful approximation. We checked that this is also the case for the POD method. Table 11 shows the computational data obtained for the reductions. We can clearly see that the number of equations of the reduced models are higher than the ones obtained for the behaviors studied in the previous, as well no speed-up is obtained and the errors are not assumable. Fig. 15 shows a comparison of the solutions obtained for the methods presented above. It is clear that the reduced solutions do not provide a good approximation of the behavior studied, and that the reduced behavior is very erratic. This loss in accuracy could be due to the fact that the original system (13) exhibits chaotic-like behavior for this range of parameters, where each cell is independent of the others and its behavior cannot be expressed as a linear combination of the behaviors of the other cells for the whole length of the integrating range. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 19 A. Bandera et al. Table 11 Comparison of the results obtained, considering system (13) for 𝑁= 100 with 𝑐𝛼= −0.5and 𝑐𝛽∈ [−0.25,0.25]. POD POBD3 POBD6 APOBD Number of equations 109 73+33+31 37+36+19+16+17+15 109 Sparsity index 0 0.6797 0.7676 0 Speed-ups <1 <1 <1 <1 Relative error 0.96 0.93 0.90 0.96 Fig. 15. Representation of the solutions obtained for the original system and the reduced ones. (a) Original, (b) POD, (c) POBD with 3 blocks, (d) POBD with 6 blocks and (e) APOBD. A.2. Limitations in the case of nonlinear coupling Another limitation could be faced if the original model exhibit too many nonlinearities. In that case, the EIM could not be able to reproduce all the variability of the nonlinear terms, and we should need to go back to the original space to compute the nonlinearities. In the following example, we present a modified version of model (9) with nonlinear coupling: 𝑂𝑖 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 𝑥𝑖=𝜏(−𝑦𝑖+𝑓(𝑥𝑖) − 𝜙𝑓(𝑧𝑖)), 𝑦𝑖=𝜏𝜀𝑘𝑖(𝑎0𝑥𝑖+𝑎1𝑦𝑖+𝑎2+1 𝑁∕2 𝑁 ∑ 𝑗=1 𝑐𝑖𝑗 𝜎(𝑥𝑖−𝑥𝑗)), 𝑧𝑖=𝜏𝜀 (𝜙𝑟(𝑥𝑖) − 𝑧𝑖−𝑧𝑏 𝜏𝑧), (19) with the same parameters as for (9), and 𝜎(𝑥) = 2 1+exp(−2𝑥)− 1. The selection of the nonlinear coupling has been made so we keep qualitatively similar behaviors as for the original model (9). In particular, we will be focusing on the AiP behavior with 𝑐𝛽= −0.25. As done for the previous test cases, we present in Table 12 the number of equations for the reduced problem obtained with each method. Furthermore, in Fig. 16, we represent the evolution of sparsity, speed-up and error, respectively. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 20 A. Bandera et al. Table 12 Comparison of the number of equations of the reduced systems for the original system (19) in AiP behavior. Number of cells POD APOBD 100 47 19+16+13+12 500 55 21+22+13+13 1000 58 23+23+14+14 2000 60 23+24+14+14 Fig. 16. Comparisons of the results obtained for each method in the AiP behavior of model (19). The left panel shows the sparsity index, the central panel shows the speed-ups, and the right panel shows the relative error. Fig. 17. Comparisons of the solutions obtained in the (SW) behavior. (a) Original. (b) POD. (c) APOBD with 3 blocks. We can see that although that we obtain better results in terms of the sparsity index and relative error in comparison with the POD method, the speed-up is similar to the POD, but very low in comparison with the previous tests. This is due to the bottleneck of computing the nonlinear coupling in the original space. A.3. Limitations in the Hindmarsh–Rose network model In the case of the Hindmarsh–Rose Network model, we have a limitation when trying to reduce the system in the Square-Wave Bursting (SW) behavior (𝐼= 3,𝑏= 2.6,𝜀= 0.01), see Fig. 17(a). In this type of bursters, the stable equilibrium corresponding to the resting state disappears through a fold bifurcation and the limit cycle attractor corresponding to the spiking state disappears through a Homoclinic bifurcation. The presence of this Homoclinic bifurcation makes the system to be near to chaotic behavior. In this case, the reduced models are not able to truly reproduce the behavior of the original system (14) as they exhibit some erratic behavior after a certain amount of time, see Figs. 17(b) and 17(c), due to the presence of the Homoclinic bifurcation. References [1] Ton R, Deco G, Daffertshofer A. Structure-function discrepancy: inhomogeneity and delays in synchronized neural networks. PLoS Comput Biol 2014;10(7):e1003736. [2] Daffertshofer A, Pietras B. Phase synchronization in neural systems. Synergetics 2020;221–33. [3] Coombes S. Next generation neural population models. Front Appl Math Stat 2023;9:1128224. [4] Touboul J. Mean-field equations for stochastic firing-rate neural fields with delays: Derivation and noise-induced transitions. Physica D 2012;241(15):1223–44. Communications in Nonlinear Science and Numerical Simulation 131 (2024) 107844 21 A. Bandera et al. [5] Faugeras O, Grimbert F, Slotine J-J. Absolute stability and complete synchronization in a class of neural fields models. SIAM J Appl Math 2008;69(1):205–50. [6] Byrne Á, Ross J, Nicks R, Coombes S. Mean-field models for EEG/MEG: from oscillations to waves. Brain Topogr. 2022;35(1):36–53. [7] Traub RD, Kopell N, Bibbig A, Buhl EH, LeBeau FE, Whittington MA. Gap junctions between interneuron dendrites can enhance synchrony of gamma oscillations in distributed networks. J Neurosci 2001;21(23):9478–86. [8] Semenov DM, Fradkov AL. Adaptive synchronization in the complex heterogeneous networks of Hindmarsh–Rose neurons. Chaos Solitons Fractals 2021;150:111170. [9] Somers D, Kopell N. Waves and synchrony in networks of oscillators of relaxation and non-relaxation type. Physica D 1995;89(1–2):169–83. [10] Drover J, Rubin J, Su J, Ermentrout B. Analysis of a canard mechanism by which excitatory synaptic coupling can synchronize neurons at low firing frequencies. SIAM J Appl Math 2004;65(1):69–92. [11] Ermentrout B, Wechselberger M. Canards, clusters, and synchronization in a weakly coupled interneuron model. SIAM J Appl Dyn Syst 2009;8(1):253–78. [12] Köksal Ersöz E, Desroches M, Krupa M. Synchronization of weakly coupled canard oscillators. Physica D 2017;349:46–61. [13] Rotstein HG, Kuske R. Localized and asynchronous patterns via canards in coupled calcium oscillators. Physica D 2006;215(1):46–61. [14] Zhabotinsky AM, Rotstein HG, Epstein IR, Kopell N. A canard mechanism for localization in systems of globally coupled oscillators. SIAM J Appl Math 2003;63(6):1998–2019. [15] Best J, Borisyuk A, Rubin J, Terman D, Wechselberger M. The dynamic range of bursting in a model respiratory pacemaker network. SIAM J Appl Dyn Syst 2005;4(4):1107–39. [16] Izhikevich EM. Synchronization of elliptic bursters. SIAM Rev 2001;43(2):315–44 [17] Coombes S, Wedgwood KC. Weakly coupled oscillator networks. In: Neurodynamics. Springer; 2023, p. 227–91. [18] Chapelle D, Gariah A, Sainte-Marie J. Galerkin approximation with proper orthogonal decomposition: new error estimates and illustrative examples. ESAIM Math Model Numer Anal 2012;46(4):731–57. [19] Freund RW. Krylov-subspace methods for reduced-order modeling in circuit simulation. J Comput Appl Math 2000;123(1–2):395–421. [20] Freund RW. Model reduction methods based on krylov subspaces. Acta Numer 2003;12:267–319. [21] Freund RW. The SPRIM algorithm for structure-preserving order reduction of general RCL circuits. Model Reduct Circuit Simul 2011;25–52. [22] Kahlbacher M, Volkwein S. Galerkin proper orthogonal decomposition methods for parameter dependent elliptic systems. Discuss Math Differ Incl Control Optim 2007;27(1):95–117. [23] Kerschen G, Golinval J-c, Vakakis AF, Bergman LA. The method of proper orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview. Nonlinear Dyn 2005;41:147–69. [24] Kostova-Vassilevska T, Oxberry GM. Model reduction of dynamical systems by proper orthogonal decomposition: Error bounds and comparison of methods using snapshots from the solution and the time derivatives. J Comput Appl Math 2018;330:553–73. [25] Kunisch K, Volkwein S. Galerkin proper orthogonal decomposition methods for parabolic problems. Numer Math 2001;90:117–48. [26] Rathinam M, Petzold LR. A new look at proper orthogonal decomposition. SIAM J Numer Anal 2003;41(5):1893–925. [27] Volkwein S. Proper orthogonal decomposition: Theory and reduced-order modelling. Lecture notes, vol. 4, (no. 4):University of Konstanz; 2013, p. 1–29. [28] Xu K-L, Jiang Y-L, Li Z, Li L. Model reduction of discrete time-delay systems based on Charlier polynomials and high-order krylov subspaces. Linear Algebra Appl 2023. [29] Benner P, Gugercin S, Willcox K. A survey of projection-based model reduction methods for parametric dynamical systems. SIAM Rev 2015;57(4):483–531. [30] Kumar R, Ezhilarasi D. A state-of-the-art survey of model order reduction techniques for large-scale coupled dynamical systems. Int J Dyn Control 2023;900–16. [31] Bandera A, Fernández-García S, Gómez-Mármol M, Vidal A. A multiple timescale network model of intracellular calcium concentrations in coupled neurons: Insights from ROM simulations. Math Model Nat Phenom 2022;17:11. [32] Anderson S, White C, Farhat C. Space-local reduced-order bases for accelerating reduced-order models through sparsity. Internat J Numer Methods Engrg 2022. [33] Barrio R, Shilnikov A. Parameter-sweeping techniques for temporal dynamics of neuronal systems: case study of Hindmarsh-Rose model. J Math Neurosci 2011;1(1):1–22. [34] Bonaventura L, Fernández-García S, Gómez-Mármol M. Efficient implicit solvers for models of neuronal networks. 2022, arXiv preprint arXiv:2210.01697. [35] Etémé AS, Tabi CB, Mohamadou A. Long-range patterns in Hindmarsh-Rose networks. Commun Nonlinear Sci Numer Simul 2017;43:211–9. [36] Hindmarsh JL, Rose R. A model of neuronal bursting using three coupled first order differential equations. Proc R Soc Lond Ser B Biol Sci 1984;221(1222):87–102. [37] De Vries G. Multiple bifurcations in a polynomial model of bursting oscillations. J Nonlinear Sci 1998;8:281–316. [38] Izhikevich EM. Neural excitability, spiking and bursting. Int J Bifurcation Chaos 2000;10(06):1171–266. [39] Pernarowski M. Fast subsystem bifurcations in a slowly varying liénard system exhibiting bursting. SIAM J Appl Math 1994;54(3):814–32. [40] Smolen P, Terman D, Rinzel J. Properties of a bursting model with two slow inhibitory variables. SIAM J Appl Math 1993;53(3):861–92. [41] Quarteroni A, Manzoni A, Negri F. Reduced basis methods for partial differential equations. Springer; 2016. [42] Fernández-García S, Vidal A. Symmetric coupling of multiple timescale systems with mixed-mode oscillations and synchronization. Physica D 2020;401:132129. [43] Pernarowski M, Miura R, Kevorkian J. Perturbation techniques for models of bursting electrical activity in pancreatic 𝛽-cells. SIAM J Appl Math 1992;52(6):1627–50. [44] Fallani FDV, Corazzol M, Sternberg JR, Wyart C, Chavez M. Hierarchy of neural organization in the embryonic spinal cord: Granger-causality graph analysis of in vivo calcium imaging data. IEEE Trans Neural Syst Rehabil Eng 2015;23(3):333–41. [45] Barrault M, Maday Y, Nguyen NC, Patera AT. An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations. C R Math 2004;339(9):667–72. [46] Grepl MA, Maday Y, Nguyen NC, Patera AT. Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations. ESAIM Math Model Numer Anal 2007;41(3):575–605. [47] Shilnikov A, Kolomiets M. Methods of the qualitative theory for the Hindmarsh–Rose model: A case study–a tutorial. Int J Bifur Chaos 2008;18(08):2141–68. [48] Bertram R, Butte MJ, Kiemel T, Sherman A. Topological and phenomenological classification of bursting oscillations. Bull Math Biol 1995;57(3):413–39. [49] Barrio R, Ibáñez S, Pérez L. Hindmarsh–Rose model: Close and far to the singular limit. Phys Lett A 2017;381(6):597–603. [50] Holden L, Erneux T. Slow passage through a hopf bifurcation: from oscillatory to steady state solutions. SIAM J Appl Math 1993;53(4):1045–58. [51] Holden L, Erneux T. Understanding bursting oscillations as periodic slow passages through bifurcation and limit points. J Math Biol 1993;31:351–65.