Simulation, Modeling, and Analysis of Soliton Waves Interaction and Propagation in CNN Transmission Lines for Innovative Data Communication and Processing

We present an innovative approach to study the interaction between oblique solitons, using nonlinear transmission lines, based on Cellular Neural Network (CNN) paradigm. A single transmission line consists of a 1D array of cells that interact with neighboring cells, through both linear and nonlinear connections. Each cell is controlled by a nonlinear Ordinary Differential Equation, in particular the Korteweg de Vries equation, which defines the cell status and behavior. Two typologies of CNN transmission lines are modelled: crisscross and ring lines. In order to solve KdV equations two different methods are used: 4th-order Runge-Kutta and Forward Euler methods.This is done to evaluate their accuracy and stability with the purpose of implementing CNN transmission lines on embedded systems such as FPGA andmicrocontrollers. Simulation/analysis Graphic User Interface platforms are designed to conduct numerical simulations and to display elaboration results. From this analysis it is possible both to identify the presence and the propagation of soliton waves on the transmission lines and to highlight the interaction between solitons and rich nonlinear dynamics. With this approach it is possible to simulate and develop the transmission and processing of information within large brain networks and high density sensor systems.


Introduction
In the scientific literature in the last 40 years, the presence of solitons in Nonlinear Transmission Lines (NLTL) has been extensively studied, both theoretically and experimentally [1][2][3][4].In the continuous case, the Korteweg de Vries (KdV) equation [5] can be easily deduced, by using the reductive perturbation methods [6].It is well known that this equation has solutions of soliton type, the behavior of which has been extensively studied both numerically [7,8] and analytically [9,10].Also some interesting related negative-order integrable equations, including the negative-order KdV equation and some associated properties, were obtained, such as the results in [11][12][13].Since the early work in NTLs [1,14,15], these solitons are still subject to studies [16][17][18][19][20]. Recently, the first robust electrical oscillator has been built, the dynamics of which is based on soliton waves [21,22].
Recent works present a completely different approach [23][24][25][26], based on a Cellular Neural Network (CNN).Introduced in 1988 by Chua and Yang [27,28], CNN consists of a lattice of individual cells, which interact with neighboring cells, usually to perform very advanced sensory recognition tasks.In general, the state of the CNN cells evolves in compliance with nonlinear Ordinary Differential Equations (ODEs), and it is well known that the behavior of such systems is very rich with striking manifestations of chaos.Nonlinear ODEs describe a set of physical and biological phenomena, fundamental in contemporary scientific research and dynamical systems theory.
The nonlinear dynamics of the CNN is extremely complex and is used in contexts of real-time image processing [29], to study the emergent behaviors of complex systems (e.g., Turing patterns), or as metamodel to simulate other mathematical systems, such as Cellular Automata (CA) [30,31].
i − 2 i − 1 i + 1 i + 2 i Figure 1: Each cell interacts in both linear (blue lines) and nonlinear fashion (red lines) with their adjacent cells but just in a linear way with neighboring not adjacent cells.
Figure 2: CNN circuits.The linear connections are made using active and passive resistances; whereas, non-linear connections are modeled using non-linear piece-wise functions.
CNNs can be easily implemented in hardware and, therefore, capable of processing large amounts of information very quickly and in parallel.They can potentially be used to build artificial organs, such as retinas or other sensory systems, to be embodied in artificial robots, to analyze vast amounts of data automatically coming from Magnetic Resonance Imaging (MRI), for patients with degenerative diseases [32].
In [33] CNN NLTL have been presented and implemented on FPGA.This FPGA system, called DCMARK (Distributed Computing Micro-ARCHitecture), allows solving complex differential equations using less time than other systems.DCMARK uses the 1D KdV equation just as a benchmark to evaluate and test the proper working of the system.Ultimately, it can produce a complex nonlinear dynamics and show soliton waves.One of the main problems in the soliton theory is the interaction between solitons in multiple dimensions [34][35][36][37], due to its connection with application.For example, this problem is of fundamental importance in optical computing [38].
This approach allows the introduction of crisscrossed transmission lines for studying the interaction among soliton waves.Unlike traditional models based on PDEs that do not permit studying the intersection of oblique solitons, this approach allows crossing the lines of propagation in a very simple way and observing the behavior of solitons, when they cross another line or when they collide with each other.Throughout this paper, we will see how this nonlinear dynamics is extremely rich and varied, as it presents unexpected and unpredictable behaviour.Such dynamics, never observed in these contexts, may be useful for engineering applications, specifically dedicated to physical systems for transmitting information.In this paper, we deal with the following: the general concept of CNN transmission and typologies of NLTL investigated in Section 2, the simulation/analysis software environments used in Section 3, the evaluation of simulation results in Section 4, and the conclusions in Section 5.

CNN Transmission Lines
As already introduced before, a CNN transmission line can be seen as a 1D cellular structure of  elements, called "cells, " which interact through both linear and nonlinear connection with their adjacent cells and through only linear connection with the other neighbours as in Figure 1.
Let us suppose that the dynamics of each cell can be modeled through the following ODE: where is a polynomial function of  and () =  is a linear function of .If we restrict our attention to polynomials of second degree, (1) can be written as where  + (with  = −2; −1; +1; +2),  + (with  = −1; +1) are coefficients.In (3),   =   = 0. Assuming that  + and  + do not depend on the th cell, that is, we consider homogeneous networks, then (3) depends on 6 parameters  −2 ,  −1 ,  1 ,  2 ,  −1 , and  1 , representing the genome of the network.Equation (3) differs from the standard CNN, because of the presence of a reaction-diffusion term and because it is a statecontrolled network.These networks can be easily implemented in hardware Figure 2, using single-cell capacitors and resistors for active and passive connections, respectively.Nonlinear connections can be modeled through Piecewise Linear (PWL) functions.
Considering parameter that depends on the number of cells and on the spatial interval, (3) becomes This is the spatial discretized version of 1D KdV equation.In fact, (4) in the continuum limit ℎ → 0 approximates to the following 1D KdV equation: Equation ( 4) is getting solving of Kirchhoff 's law as in [25].
Simulating (4) and assuming an initial function such as a Gaussian function, we observe the emergence of traveling waves as shown in Figure 3.
Simulating a CNN transmission line of 300 cells and fixing properly the boundary conditions, traveling soliton waves are shown.In fact, as the classical solitons of the KdV equation, they have very well-defined profiles, with a particular relationship between amplitude, width, and velocity of propagation.The taller and slimmer solitons travel faster than the lower and wider ones, and they maintain their distinctive character in the interaction with the others, as it will be explained later.
If we consider   < 0 in (4), a new linear damping term is added.In this case, this equation becomes a modified KdV equation, due to the presence of inhomogeneity.If   = 1/2, (4) becomes, in the continuous case, just the cylindrical KdV equation.
Model (1), whose corresponding patterns are plotted in Figure 3, gives rise to a rich environment, in which the nonlinear dynamics can be studied.
On the basis of the same reasoning, with some modifications, it is possible to simulate the crisscross of two 1D CNN transmission lines analysing the interaction between solitons belonging to different lines.Even if the interaction of oblique solitons is a classic problem, it is still far from being solved.The main obstacle to its solution is that the equations at its basis are rare, and no approximation seems to be convincing.The interest of researchers in this area is very high, especially in the context of optical computing, because these studies may provide a significant input for new forms of computation.Now, numerical models based on CNN transmission lines can provide new and natural approaches to solve this long-standing problem.
Considering two transmission lines that intersect in a given cell (, ), called cross cell, as in Figure 4, the state of each cell, except for cross cell (, ), is governed by (4).Not only does the neighborhood of cross cell (, ) consist of the two cells at its left and right, respectively but also it is made up of its upper and lower cells.As before, the neighboring cells, up and down, have a linear excitatory synaptic connection, whereas the other connection is nonlinear.Upper or lower cells have rather inhibitory synaptic connections.
From a circuital point of view, this leads to a configuration similar to that in Figure 5.
All the cells of the two CNN transmission lines have (4) as spatial discretized state equation; just at the point of intersection between the two lines, cross cell (, ) is ruled by the following spatial discretized state equation: This equation that governs the behavior of the cross cell could be obtained as done for (4) in [25], solving Kirchhoff 's law for the circuit in Figure 5 and finding the expression of current through the capacitor (, ).The state (, ) of the cross cell and the voltage across the capacitor is represented by , the inhibitory synaptic connection,   , the excitatory synaptic connection, and   , the coefficient of the nonlinear term.Then the equation of the state of the cross cell (, ) becomes from which taking   = /2,  = 2ℎ 2 , and    = 2ℎ/3, where ℎ = cost, and grouping together equation terms state equation ( 6) is obtained.ℎ is a normalization parameter by which it is possible to tune the system model to improve convergence and stability.The dynamics of cross cell is regulated by the state of eight neighboring cells, with radius 2. Neighboring cells, with radius 1, interact with the state of the cell in both linear and nonlinear ways.
Hence the soliton interaction is studied in two different scenarios: a ring line on which there are a single soliton propagating and a crisscross of two lines with a soliton on both lines.The motivation for the choice of these two topological setups sprang from facility for analysis and detection of soliton dynamics.In particular, the second one (relating to crisscross setup) is the best for analyzing the interaction between oblique solitons, very useful for highlighting typical dynamical behaviors of interaction, without any complications related to topology.
In order to simulate the presence of a soliton, we can use, as an initial condition, a function such as a Gaussian function or a square hyperbolic secant function with different magnitudes.
Considering a square hyperbolic secant function as initial condition: This function avoids divergence integration problems, thanks to its zero-tangent envelope for  → ±∞.
According to [7] if one considers the following initial condition: the eigenvalues of Schrodinger equation spectrum to which every Korteweg de Vries solution is related are So, if  is an integer number, there will be the birth of  solitons from the initial state, otherwise there will be  solitons with a radiation tail coming after.
The magnitude of solitons is and the velocity of propagation is As already introduced, amplitude, width, and velocity are connected.For example, amplitude is directly proportional to velocity of propagation.Higher amplitude corresponds to higher velocity of propagation.
In order to increase slightly the stability and accuracy of Forward Euler method which will be firstly validated using a custom MATLAB tool and then implemented on embedded platforms, a particular handling is done on discretized state equations.
Starting from (4) and sorting the equation terms,  single equations [39] are obtained from the KdV equation: where  = 0, . . .,  is the space iteration index and Δ is the spatial step of the discrete grid.This numerical discretization of spatial derivative terms of ( 5) has been done using a spacecentered finite difference method [40].
For the time derivative term of ( 5), just for the first iteration, we used a forward-time finite difference method as in [7,41] because there is no preceding value at the first step of numerical integration process, obtaining (13).Hence, for the other iterations, we used a centered-time finite difference method, obtaining (14).We set where  = 0, . . .,  are the time iteration index,  = 0, . . .,  are the spatial iteration index, and Δ is the integration time.
The same considerations about the spatial and time discretization methods have also been done for the crisscrossed lines scenario as well as the ring line scenario.Starting from (6) and sorting the equation terms, the spatial discretized state equation for the cross cell of intersection between the two lines is obtained: where ,  = 0, . . .,  are the spatial iteration index for the two lines and Δ is the space step of the discrete grid.Equation ( 16) is the equation in which the forward-time finite difference method is applied while in ( 17) the centred-time difference method is used, for time discretization: where  = 0, . . .,  is the time iteration index, ,  = 0, . . .,  are the spatial iteration index for the two lines, and Δ is the integration time.
Using this combined approach a stable propagation of a soliton through all cells for all time cycles is obtained.This kind of discretization is less accurate than other types, but it is also the best technique in terms of implementation easiness and resources saving on embedded systems.

Simulation/Analysis Settings and Environment
In the literature, there are many numerical methods to solve state ODEs.In this work, two methods are used for two different aims: a 4th-order Runge-Kutta (RK4) method in order to conduct accurate high level analysis and a simpler Forward Euler (FE) method to be implemented, after a validation phase, on embedded platforms such as FPGAs, microcontrollers, DSP, or ASIC.FE method is chosen in order to relax the computing load of embedded platforms to the detriment of accuracy and stability.
With the goal of conducting extensive simulations and analysis, two Graphic User Interface (GUI) environments have been designed (SIMulation-ANalysis-SOLiton tools).The first, called C++ SIMANSOL, has been designed using C++ language (Figure 6); the latter, called ML SIMANSOL [42], instead has been designed in MATLAB environment (Figure 7).C++ SIMANSOL tool has been used mainly for fast and accurate analysis while ML-SIMANSOL tool has been used for validation and investigation in order to verify the efficiency, accuracy, and stability of FE method.MATLAB is preferred to other high level languages because of its optimal aptitude to handle matrices and vectors easily.Having these two high level environments allows both to compare their respective simulation results and to validate the ODE solving process.A comparison between Forward Euler method and Runge-Kutta method is not to show an evident difference between these two methods trivially but just for verifying if the stability and accuracy of Forward Euler method are acceptable using a limited number of CNN cells with respect to Runge-Kutta method.The impossibility of implementing on FPGA one 1D CNN, formed by a huge number of computing processors (one for each cell), causes problems of equation solution instability and divergence.So, it was very important to understand the minimum number of processors to implement on FPGA in order to minimize these problems.In addition, implementing a Forward Euler method for numerical integration permits saving a lot of FPGA resources guaranteeing the possibility to implement larger and larger Cellular Neural Network on the FPGA [26], allowing facing more complex dynamics problems.
In the C++ SIMANSOL tool (Figure 6), there are three different windows: initial input function window, simulation setup window, and multiwires grid setup window.It is possible to control all simulation parameters such as integration method (Euler, RK4th, RK2th, RK45, and RKAdaptive), boundaries (periodic, zero flux, and fixed), number of cells, number of iterations, initial input function (Gaussian, hyperbolic secant, sin, soliton, impulse, etc.), time step, spatial step, and simulation scenario (single line, crisscrossed lines, and grid network).
Instead MATLAB SIMANSOL tool (Figure 7) consists of only one main window from which it is possible to control several simulation parameters such as number of cells, number of iterations, soliton amplitude, time step, spatial step, analysis type (time, spatial or time/spatial), simulation scenario (ring line or crisscross lines), and plot simulation results.This tool, as already said, implements the FE integration method.
Using these tools it is possible to give a complete vision and to understand complex dynamics phenomena.Three types of graphs can be displayed: a state cell variable versus time graph, in which the time evolution of the state variable associated with a certain cell of the ring line is shown; a state cell variable versus space graph, where the value of all state cell variable associated with all cells is displayed for a well-defined time step; a time/space graph, in which the time and space evolution of the state cell variable of all cells is highlighted.This last graph is very interesting because it allows getting a global vision about the propagation of solitons and their interaction.

Ring Line Analysis.
In the first scenario, a ring line of 100-400 cells is considered according to what was said previously.The ring-fashion feature is obtained simply by imposed periodic boundary conditions.Three MATLAB simulations using the MATLAB (ML) SIMANSOL tool have been conducted under the following conditions (Table 1).

1st Test: Soliton with Magnitude 2.
In Figure 8 the state variable of 10th cell in function of time (number of iterations) is shown.A single soliton with magnitude 2 propagates  through the ring line performing several loops.This soliton has a velocity of propagation of 4 according to (9), completing a single loop in about 1250 iterations or rather 12.5 s.There are not any numerical instability or divergence problems during simulations.This simulation allows verifying the quality of Forward Euler method to be used for implementation of soliton ring line on embedded devices [26,33].As shown in Figure 9, the narrow parallel lines with slope of 4 represent the single soliton which performs loops through the ring.
As already explained in Section 2, with an amplitude of 2, according to [7], just a single soliton is foreseen.

2nd
Test: Soliton with Magnitude 6.This time the state cell variable of 200th cell is analysed (Figure 10).The soliton now having an amplitude of 6 generates two solitons.Hence at the 150000th iteration, it is possible to see the interaction between the two solitons with magnitudes 8 and 2, respectively.These two solitons come from a single soliton with initial magnitude 6 which propagates through the ring line.They have a velocity of propagation of 16 and 4, respectively.This interaction is possible because of different velocity of two solitons.
From Figure 11 it is possible to understand the velocity of solitons on the basis of straight lines.The narrow parallel lines with slope of 16 represent the faster soliton with magnitude of 8 while the large line with slope of 4 represents the soliton with magnitude of 2. In this graph, the soliton interaction  points are evident.When there is an interaction between solitons, it is possible to see a phase shift that is shown in the graph as a shift of lines.The slow soliton has a phase shift larger than fast soliton, and the direction of phase shift is positive for slow soliton.It is also interesting to see in Figure 12 the generation and propagation through the ring line of the two solitons coming   from a single soliton (blue) with magnitude 6.It is possible to see the soliton generation at five different time iterations.

3th
Test: Soliton with Magnitude 12.In this case, following the time evolution of state cell variable of 150th cell more interactions between solitons can be seen (Figure 13).At the 110000th iteration, there is an interaction between the soliton with magnitude 18 and the soliton with magnitude 2; at the 150000th iteration, it is possible to see the interaction between the two solitons with magnitudes 8 and 2, respectively, while at the 270000th iteration we find an interaction between solitons with magnitudes 18 and 8, respectively.These three solitons come from a single soliton with initial magnitude 12 which propagates through the ring line.These solitons have a velocity of propagation of 36, 16, and 4, respectively.
Also in this case from Figure 14 it is interesting to see the different velocity of solitons and the soliton interaction points.The narrow parallel lines with slope of 36 represent the faster soliton with magnitude of 18 while the other lines with slope of 16 and 4 represent the solitons with magnitude of 2 and 8, respectively.In this graph, the soliton interaction points are evident.When there is an interaction between solitons, it is possible to see a phase shift that is shown in the graph as a shift of lines as in the other cases seen previously.
From Figure 15 the generation and propagation through the ring line of three solitons coming from a single soliton  (blue) with magnitude 12 are shown.It is possible to see the soliton evolution at five different time iterations as in the previous cases.

Crisscross Lines Analysis.
In this section some simulations of crisscross transmission lines are shown.Also for these simulations the SIMANSOL tools are used.In both tools periodic boundary conditions are chosen, building a sort of crossed rings.Three different types of test are conducted.The initial conditions for these simulations are summarized in Table 2.
With the C++ SIMANSOL tool, solitons that propagate along the lines before and after the crossing point are shown on the graph while using ML-SIMANSOL tool the Space/ Time graphs for lines  and  are plotted.

1st
Test: One Soliton on Both Lines  and .As it can be seen in Figures 16 and 17, two solitons belonging to different lines meet in the crossing point and separate again, propagating on their own line as if this interaction never happened.After the interaction, they preserve their distinctiveness.The only observed phenomenon is the presence of a small dispersive track and reflected solitons, linked to the periodic boundary conditions, rather than to the interaction between the two solitons.

2nd
Test: Just One Soliton on Line .The second experiment involves the presence of a single soliton passing through the points of intersection between the two lines.When the soliton reaches the point of intersection, it slows down, while on the two lines disturbances occur.After crossing the interaction point, two new solitons emerge, one for each line; it is possible to see two wave trains as well, one for each line, but traveling in the opposite direction compared to solitons, as shown in Figures 18 and 19.
There is also a decreasing of magnitude of the two solitons; on line , after the crossing point, soliton has magnitudes 0.6 on line  and 0.5 on line  for simulation with C++ SIMANSOL tool, while 1.3 and 0.5, respectively, on average for simulation with ML-SIMANSOL tool.The larger soliton is the one that resides on the same line as the one of the incoming soliton.In addition, it is possible to see phenomena of dispersive tracks and reflected solitons too.

3th
Test: One Soliton on Line  and One Delayed Soliton on Line .The third experiment concerns the interaction between two solitons of the same amplitude but phaseshifted.The first soliton reaches the intersection cell before the other one and undergoes a process of division, similar to that described in the previous experiment.When the second soliton reaches the point of intersection, it moves in an unstable environment.The second soliton undergoes a process of division and, at the end of the interaction, a total of 4 new solitons emerge.Figures 20 and 21 show the second soliton while going through the intersection cell.Again, we note that there is no symmetry between the two pairs of solitons.On the line of the soliton, which first crosses the intersection cell, an emerging pair is generated, which has a width slightly larger than the second pair.
The generation of two small solitons and dispersive tails on both lines is also evident.
The results that emerge from this series of experiments are extremely interesting.Using both simulation tools, it is possible to see the same results, even if there are numerical mismatches because of the two different integration methods.They show that the interaction of oblique solitons is very rich and diverse.In Table 3, the results of the experiments are shown.
In the first case, two identical solitons interact and emerge unchanged after crossing.In the second case, a soliton is split into two smaller and different solitons.In the third case, two identical but phase-shifted solitons splits into four solitons, slightly different from each other.This kind of results, even if it is still under research, could give, also, some ideas about a possible utilization of these CNN transmission lines in terms of digital logic gates, making the most of their nonlinear and innovative features.

Conclusions
In this work, a new approach to the transmission lines based on nonlinear CNN is presented.After providing the basic equation that governs the status of each cell, the presence and the interaction of solitary waves have been detected.These waves interact with each other, while maintaining their distinctiveness.This behavior is not surprising, because in the case where the network is continuous the equation of state just becomes the well-known KdV equation.The introduction of nonhomogeneous terms provides a new approach to the study of wave propagation in multidimensional nonlinear and dispersive media.Furthermore, this approach allows dealing very effectively with the problems related to interaction between oblique solitons.Several simulations in two different scenarios (ring line and crisscross lines) have been conducted with the help of two custom simulation tools.The possibility to implement a simple Forward Euler method for solving KdV ODE on embedded systems has been validated.The results have great potential applications, especially in the construction of information transmission lines, embodied into other artificial intelligent media, such as automatic expert systems and robotic agents.There are many advantages using CNNs for designing transmission lines.CNNs can be easily implemented in silicon performing multiple processes in parallel, at very high speed.They constitute a very flexible simulation environment by which we investigate nonlinear dynamics and carry out a huge number of trials with different initial data.Using traditional systems, the direct observation of solitons behavior could be very complicated.Such an approach may guarantee a technological environment for developing new and unconventional forms of computation, which exploit fully nonlinear dynamics, opening interesting perspectives for the study of multidimensional solitons, with reference to optical computing applications.In addition, this innovative approach allows simulating neuromorphic systems, creating new simulation hardware for the development of brain-like physical agents.

Figure 3 :
Figure 3: The initial Gaussian data decomposed into a series of localized travelling waves.

Figure 4 : 2 Figure 5 :
Figure 4: Block diagram representation of the crisscrossed lines structure at the point of intersection between the two lines (red and yellow lines for nonlinear interactions, blue and violet lines for linear interactions).

Figure 7 :
Figure 7: MATLAB SIMANSOL GUI tool for simulation and analysis.

Figure 8 :
Figure 8: A single soliton with magnitude 2 which propagates through the ring line in correspondence with the 10th cell.

Figure 9 :
Figure 9: A Space/Time graph in which a single soliton with magnitude 2 propagates through the ring line.

Figure 10 :
Figure 10: The two solitons which propagate and interact through the ring line in correspondence with the 200th cell.

Figure 11 :Figure 12 :
Figure 11: A Space/Time graph in which a single soliton with magnitude 6 propagates through the ring line and is divided into two solitons.

Figure 13 :
Figure 13: The three solitons which propagate and interact through the ring line in correspondence with the 150th cell.

Figure 14 :
Figure 14: A Space/Time graph in which a single soliton with magnitude 12 propagates through the ring line and is divided into three solitons.

Figure 15 :
Figure 15: The generation and propagation through the ring line of three solitons from a single soliton with magnitude 12.

Figure 16 :Figure 17 :Figure 18 :
Figure 16: Simulation with C++ SIMANSOL tool: two identical solitons meet and cross with no interaction.

Figure 19 :Figure 20 :
Figure 19: Simulation with ML-SIMANSOL tool: Space/Time graphs for line  (left) and line  (right).Only one soliton on line .

= 1 )Figure 21 :
Figure 21: Simulation with ML-SIMANSOL tool: Space/Time graphs for line  (left) and line  (right).One soliton on line  and one soliton with a negative spatial delay of 25 cells on line .

Table 1 :
Ring line simulation settings.