I work in the general areas of nonlinear dynamics, time series analysis and applied stochastic dynamics. I apply my work in areas like neuroscience, geophysics and accelerator physics.
A more detailed description of my research work is given below.
Extracting information flow in networks of coupled dynamical systems from the time series measurements of their activity is of great interest in physical, biological and social sciences. Such knowledge holds the key to the understanding of phenomena ranging from turbulent fluids to interacting genes and proteins to networks of neural ensembles. Granger causality has emerged in recent years as a leading statistical technique for accomplishing this goal. The definition of Granger causality is based on the theory of linear prediction and its original estimation framework requires autoregressive (AR) modeling of time series data. This is a parametric approach to estimating Granger causality.
The above parametric approach has certain weaknesses. We overcame this by proposing a nonparametric approach to estimate Granger causality (in the spectral domain) directly from Fourier and wavelet transforms of data, eliminating the need of explicit AR modeling. This approach makes use of spectral factorization. Time-domain Granger causality can be obtained by integrating the corresponding spectral representation over frequency. This method was successfully applied to both simulated data and data from a neuroscience experiment.
Simultaneous recordings of spike trains from multiple single neurons are becoming commonplace. Mathematically, spike trains can be modeled as point processes. Understanding the interaction patterns among these spike trains (point processes) remains a key research area. A question of interest is the evaluation of information flow between neurons through the analysis of whether one spike train exerts causal influence on another. For continuous-valued time series data, Granger causality has proven an effective method for this purpose. However, the basis for Granger causality estimation is autoregressive data modeling, which is not directly applicable to spike trains. Various filtering options distort the properties of spike trains as point processes. We proposed a new nonparametric approach to estimate Granger causality directly from the Fourier transforms of spike train data. We validated the method on synthetic spike trains generated by model networks of neurons with known connectivity patterns and then applied it to neurons simultaneously recorded from the thalamus and the primary somatosensory cortex of a squirrel monkey undergoing tactile stimulation.
Computing Granger causal relations among bivariate experimentally observed time series has received increasing attention over the past few years. Such causal relations, if correctly estimated, can yield significant insights into the dynamical organization of the system being investigated. Since experimental measurements are inevitably contaminated by noise, it is thus important to understand the effects of such noise on Granger causality estimation. Using the theory of stochastic processes, we showed mathematically that two effects can arise: (1) spurious causality between two measured variables can arise and (2) true causality can be suppressed. We also demonstrated these effects numerically. We then provided a denoising strategy to mitigate this problem. Specifically, we proposed a denoising algorithm based on the combined use of the Kalman filter theory and the Expectation-Maximization (EM) algorithm. Numerical examples were used to demonstrate the effectiveness of the denoising approach.
We applied the above method to denoise two datasets of local field potentials recorded from monkeys performing a visuomotor task. For the first dataset, it was found that the analysis of the high gamma band (60-90 Hz) neural activity in the prefrontal cortex is highly susceptible to the effect of noise, and denoising leads to markedly improved results that were physiologically interpretable. For the second dataset, Granger causality between primary motor and primary somatosensory cortices was not consistent across two monkeys and the effect of noise was suspected. After denoising, the discrepancy between the two subjects was significantly reduced.
More than 99% of the species that ever existed on the surface of the earth are now extinct and their extinction on a global scale has been a puzzle. One may think that a species under an external threat may survive in some isolated locations leading to the revival of the species. Using a general model we show that, under a common external forcing, the species with a quadratic saturation term first undergoes spatial synchronization and then extinction. The effect can be observed even when the external forcing acts only on some locations provided the dynamics contains a synchronizing term. Absence of the quadratic saturation term can help the species to avoid extinction.
Among the awesome repertoire of tasks that the human brain can accomplish, one of the more fascinating ones is how the electrical activity of millions of brain cells (neurons) is translated into precise sequences of movements.
One of the greatest challenges in applied neuroscience is to build prosthetic limbs controlled by neural signals from the brain. The ultimate goal is to provide paralytic patients and amputees with the means to move and communicate by controlling the prosthetic device using brain activity. Scientists and engineers are slowly getting closer to building such devices thanks to studies revealing a strong connection between the activity of neurons in the brain’s cerebral cortex and the movements of limbs.
With the above goal in mind, we have performed some preliminary analysis of experimental data from a macaque monkey. We considered an experiment where the monkey moves its hand in one of the eight directions based on visual cues. The goal was to predict the direction of movement correctly before the hand movement actually starts. Signals were recorded from an array of microelectrodes implanted in the monkey’s motor cortex. As a first step, we computed the summed power in the gamma frequency band for each channel and trial. Using the Mahalanobis distance, we then predicted the direction of movement from a single trial multichannel data. We obtained successful predictions significantly above the random chance level.
Presently we are conducting a major brain-machine interface experiment and results are expected soon.
Given the deluge of multi-channel data generated by experiments in neuroscience, the role of multivariate time series processing, especially the nonlinear time series processing, has become crucial. In particular, obtaining causal relations among signals is important in identifying functional relations between different regions of the brain. Over the years Gersch causality (which uses partial coherence analysis) has been employed explicitly or implicitly by many researchers as a way of identifying connectivity, sources of driving or causal influence.
For real world neurobiological data, a complicating factor is that the time series picked up by a sensor is inevitably the mixture of the signal of interest (e.g. theta oscillation in the hippocampus) and other unrelated processes, collectively referred to as the measurement noise. We studied the effectiveness of the two different causality measures, the Gersch causality and the Granger causality, in the above noisy situation.
We studied experimental recordings from the limbic system of the rat during theta oscillations and showed that the application of partial coherence and Gersch causality can lead to contradictory conclusions. We hypothesized that the observed phenomena can be explained by the differential levels of signal to noise ratio. We performed extensive simulations involving multiple time series to confirm this. The main point of our work was that, despite its wide use, partial coherence based Gersch causality is extremely sensitive to measurement noise and often could lead to erroneous results. We then applied Granger causality to the data and simulation results and showed that it is robust against moderate levels of noise.
We addressed the problem of obtaining nonlinear causality measures. We extended the Granger causality theory to nonlinear time series by incorporating the embedding reconstruction technique for multivariate time series. A three-step algorithm was presented and used to analyze various nonlinear time series. We demonstrated that our method always gives reliable results.
When three time series have to be analyzed, the conditional nonlinear causality index proposed by us can distinguish between direct and indirect causal relationships between any two of the time series. This is not possible using simple pairwise analysis.
For oscillatory neural networks, the stability of the equilibrium point is an important property for many neural computational applications. Aided by results from nonnegative matrices we demonstrated, using two simple models of coupled neural populations of arbitrary size, a general methodology that can yield explicit bounds on the individual coupling strengths to ensure the stability of the equilibrium point.
First passage time problem for stochastic processes is a well known problem with applications in physics, chemistry, biology and engineering. The exact first passage time distribution is known only for a few specific stochastic processes (the most notable being the Brownian motion). Finding the exact first passage time distribution for a broader class of stochastic processes is therefore an important problem. We obtained the exact first passage density function for Levy type stochastic processes with zero and nonzero drifts. This work involved application of results from probability theory, fractional calculus and Fox or H-functions. Further, asymptotic results for another class of stochastic processes (fractional Brownian motion and its generalizations) were also obtained.
Coupled dynamical systems are used to model spatially extended nonlinear systems (such as the brain) by utilizing well understood dynamical systems as building blocks. In coupled dynamical systems, the individual constituent systems (called nodes) are coupled to one another using specified coupling strengths. A very interesting form of behaviour which occurs in coupled dynamical systems is the synchronization of the nodes.
This has found applications in a variety of fields including communications, optics, neuroscience, neural networks and geophysics. An essential prerequisite for these applications is to know the bounds on the coupling strengths so that the stability of the synchronous state is ensured. Previous attempts aimed at obtaining such conditions have typically looked either at systems of very small size or at very specific coupling schemes (diffusive coupling, global all to all coupling etc. with a single coupling strength).
We obtained stability bounds for the synchronized state in coupled systems with arbitrary coupling. These bounds extend previously available results. Based on the Gershgorin disc theory and the stability region/master stability function, we developed general approaches that yield constraints directly on the coupling strengths which ensure the stability of synchronized dynamics.
In a reaction-diffusion system, at the Turing bifurcation point, the global equilibrium state becomes unstable and an inhomogeneous spatial pattern known as a Turing pattern emerges. Following Turing’s classic work, pattern formation in reaction-diffusion systems has become a major topic of investigation both theoretically and experimentally with applications in diverse fields of science and engineering.
We derived explicit analytical expressions defining the stability region in the parameter space spanned by the coupling strengths. We demonstrated that, by following different paths in the parameter space, different spatiotemporal patterns can be selectively realized. Conversely, given a desired spatiotemporal pattern, we were able to design coupling schemes which give rise to that pattern as the coupled system evolves. Although the spatial patterns are the same for a given coupling matrix, we showed that, depending on whether the original synchronized state is a fixed point or a chaotic attractor, the temporal evolution of the patterns is either constant in time or modulated in an on-off intermittent fashion.
To assess whether a given time series can be modeled by a stochastic process possessing long memory one usually applies one of two types of analysis methods: the spectral method and the random walk analysis. We showed that each one of these methods used alone can be susceptible to producing false results. We thus advocated an integrated approach which requires the use of both methods in a consistent fashion. We provided the theoretical foundation of this approach and illustrated the main ideas using examples. We also analysed the observation of long range anticorrelation (Hurst exponent $H<1/2$) in real world time series data. The very peculiar nature of such processes was emphasized in light of the stringent condition under which such processes can occur. Using examples we discussed the possible factors that could contribute to the false claim of long range anticorrelations and demonstrated the particular importance of the integrated approach in this case.
Chaotic dynamics is an important field of research with applications in many areas of science and engineering. Specific examples of such areas include astrophysical, biological and chemical systems, mechanical devices, models of the weather, lasers, plasmas and fluids, to mention a few. It is important to be able to identify chaotic systems since one can predict the future behaviour of such systems only for a limited period of time. This is due to the fact that these systems exhibit sensitive dependence on initial conditions, i.e. two trajectories which start out close to each other diverge exponentially. The Lyapunov exponents precisely quantify this exponential divergence and a positive Lyapunov exponent implies the presence of chaotic behaviour in the system. Moreover, practical algorithms to compute Lyapunov exponents can be devised. For these reasons, Lyapunov exponents have become the primary tool to detect chaos in a given system. Any good method for computing Lyapunov exponents has, therefore, wide applicability in several areas of science and engineering.
Chaotic dynamics is an important field of research with applications in many areas of science and engineering. Specific examples of such areas include astrophysical, biological and chemical systems, mechanical devices, models of the weather, lasers, plasmas and fluids, to mention a few. It is important to be able to identify chaotic systems since one can predict the future behaviour of such systems only for a limited period of time. This is due to the fact that these systems exhibit sensitive dependence on initial conditions, i.e. two trajectories which start out close to each other diverge exponentially. The Lyapunov exponents precisely quantify this exponential divergence and a positive Lyapunov exponent implies the presence of chaotic behaviour in the system. Moreover, practical algorithms to compute Lyapunov exponents can be devised. For these reasons, Lyapunov exponents have become the primary tool to detect chaos in a given system. Any good method for computing Lyapunov exponents has, therefore, wide applicability in several areas of science and engineering.
Hamiltonian systems form an important class of dynamical systems and can be used to describe many physical systems. An important question concerning any such system is its long-term nonlinear stability. Specifically, we are interested in the long-term stability of particles being transported through this system. Generically these systems are non-integrable and it is therefore very difficult to give stability criteria in an analytic form. A possible solution is to numerically follow the trajectories of particles through the system for a large number of periods (a process that goes by the name of tracking). One could then attempt to infer the stability of motion in the system by analyzing these tracking results.
The most straightforward method that can be used to perform this long-term tracking is numerical integration. However, this method is too slow for analyzing the stability of very complicated systems. Further, any method that we use should preserve the symplectic structure that is inherent in any Hamiltonian system. Conventional numerical integration methods, like the Runge-Kutta method, do not preserve this structure. This can lead to either spurious chaotic behaviour or damping. Consequently, there is a possibility of getting wrong stability results when using such non-symplectic integrators. Therefore, we need a symplectic integrator that is both fast and accurate.
Several symplectic integration methods have been discussed in the literature using generating functions. These are typically implicit methods and employing these methods requires one to use Newton’s method with its attendent questions of convergence etc. We have obtained a new method for symplectic integration of nonlinear Hamiltonian systems. In this method, the Hamiltonian system is first represented as a symplectic map. This map (in general) has an infinite Taylor series expansion. In practice, one can compute only a finite number of terms in this series. This gives rise to a truncated map approximation of the original map. This truncated map is, however, not symplectic and can lead to wrong stability results when iterated (as noted earlier). To overcome this problem, the map is refactorized as a product of special maps called “jolt maps” in such a manner that symplecticity is maintained. These jolt maps are nilpotent operators of rank 2 and hence have only two nonzero terms in their Taylor series expansion when acting on phase space variables. The method makes use of properties of the Lie group SU(3) and its associated discrete subgroups. Moreover, a scheme to optimize the number of jolt maps required in this method have been outlined. This involves splitting the integration over SU(3) into an integration over the group SO(3) followed by an integration over the coset space SU(3)/SO(3).
We have applied the above symplectic integration method using jolt factorization to the symplectic map describing the nonlinear pendulum Hamiltonian. This example was chosen since the pendulum is the prototypical nonlinear Hamiltonian system which has the further advantage that analytical solutions are available facilitating easy comparisons. The symplectic map corresponding to the pendulum Hamiltonian was obtained and truncated at the eighth order. This map was then refactorized using jolt maps. Using this for integration, results which agreed well with the exact results were obtained. In contrast, a non-symplectic eighth order integration method was shown to give rise to spurious chaotic behaviour. This demonstrated the advantages of using the symplectic integration method outlined above.
Using solvable and polynomial symplectic maps, symplectic integration algorithms which can be easily generalized to higher dimensions were also obtained.
In the study of nonlinear Hamiltonian dynamics, the real symplectic group Sp(2n,R) and its compact subgroups play an important role. Quite often, one studies the single particle dynamics of nonlinear Hamiltonian systems. Since this has three degrees of freedom, the relevant group is Sp(6,R). Moreover, the equations of motion are formulated in terms of phase space variables (generalized coordinates and momenta). In particular, in Lie perturbation theory of Hamiltonian dynamics, homogeneous polynomials of phase space variables play a central role. Therefore, it is important to study the representations of Sp(6,R) carried by these polynomials. Further, in deriving symplectic integration algorithms for Hamiltonian systems (which was discussed in the previous section), representations of compact subgroups of Sp(6,R) (especially SU(3)) carried by homogeneous polynomials are required. They may also be useful in deriving metric invariants for symplectic maps.
For the above reasons, we studied representations of Sp(6,R) and SU(3) carried by homogeneous polynomials of phase space variables. It was proved that homogeneous polynomials of degree m carry the irreducible representation (m,0,0) of both the Lie algebra sp(6,R) and the corresponding Lie group Sp(6,R). Relations between representations of sp(6,R) and Sp(6,R) were established. Under the action of SU(3), the representations obtained above are reducible and can be written as a direct sum of irreducible representations of SU(3). Finally, explicit expressions for SU(3) states within these representations were given in terms of phase space variables. Such expressions are useful in symplectic integration theory.
We also obtained a rather unconventional real basis for the symplectic algebra sp(2n,R). We demonstrated the utility of this basis for practical computations by giving a simple derivation of the second and fourth order indices of irreducible representations of sp(2n,R).
Understanding climatic dynamics and variability is an important problem in geophysics. One can attempt to solve this problem by developing elaborate climate models and studying them numerically. This requires tremendous amount of computer time and effort. Therefore, it would be good to have some simple measures of climatic variablity which can be used to understand the underlying dynamics in a gross sense and to exhibit linkages between different climatic dynamics.
We proposed a novel climate predictability index. It quantifies, in a simple manner, the predictability of three major components comprising the climate – temperature, pressure and precipitation. The quantification was done using a fractal dimensional analysis of the corresponding time series as described below. It is well known that a geophysical time series can be modelled as a fractional Brownian motion. Using this fact, it is possible to apply rescaled range (R/S) analysis to the time series and calculate its so-called Hurst exponent. This Hurst exponent is related in a simple manner to the fractal dimension of the time series curve. Fractal dimensions close to 1.5 and 1.0 correspond to low and high predictability respectively. This correspondence was exploited to define a predictability index. The climate predictability index was calculated for stations in India using the Global Historical Climatology Network dataset. Change in the index with the seasons suggested a strong influence of more than one climatic dynamics. In such cases, calculations done using mean yearly data were shown to be suspect. The index was shown to be useful in studying the interplay between various climatic components (viz. temperature, pressure and precipitation). Changes in predictability indices for temperature and pressure were seen to affect the index for precipitation. The index can be used as a discriminant for determining which stations are selected for use in developing regional climatic models.
One of the striking phenomena exhibited by a wide variety of complex adaptive systems is that individual agents or components of the system evolve to perform highly specialized tasks, and at the same time the system as a whole evolves towards a greater diversity in terms of the kinds of individual agents or components it contains or the tasks that are performed in it. Some examples of this include living systems which have evolved increasingly specialized and diverse kinds of interacting protein molecules, ecologies which develop diverse species with specialized traits, early human societies which evolve from a state where everyone shares in a small number of chores to a state with many more activities performed largely by specialists, and firms in an economic web that explore and occupy increasingly specialized and diverse niches. Since the above phenomenon is so widespread, it would be useful to have a mathematical model which captures this.
We have proposed a mathematical model of economic communities that exhibits these twin evolutionary phenomena of specialization and diversity. The community consists of N interacting agents. Each agent’s activity consists of an individual mix out of a fixed set of strategies. Agents receive a payoff depending on activities of the entire community, and independently modify their own mix of strategies to increase their payoff. The dynamics obtained is a generalized replicator dynamics and consists of a set of coupled nonlinear ordinary differential equations describing the time evolution of the activities of individual agents. We proved certain theorems and obtained numerical results for the attractors of the system. We focused our attention on attractors that have the property of individual specialization and global diversity. Under certain conditions, the desired attractors were shown to exist and have basins of attraction that cover the entire configuration space. Thus the evolutionary phenomena mentioned above occur generically in the model. We have described properties a new strategy should have (in the context of already existing activities) in order for it to be “accepted” by the community. Thus the model suggests a natural mechanism for the emergence of context dependent innovations in the community. It was also shown that under certain generic conditions on the payoff matrix parameters, the agents exhibit a collective behaviour, and for sufficiently large N, the community exhibits diversity and self-organization.
One of the classical problems in game theory and economics is that of equilibrium selection. In a typical situation, one has many candidate equilibria and the problem of selecting the right one has to be solved. One approach towards this problem is to construct dynamic models of disequilibrium wherein the desired equilibrium arises as an attractor of the proposed dynamics. The dynamics are intended to model economic agents in an interactive environment where they `learn’ from experience and adapt their strategy in real times. Exisitng models fall into two distinct classes representing two different modes of learning. In the first class of models, one retains the identity of individual agents in a finite collection and postulates adaptation rules for individual learning which lead to the aggregate behaviour being studied. In the second class of models, the population itself adapts by incrementally reinforcing or discouraging (implicitly, through birth-death mechanisms) the strategies that give respectively higher and lower payoffs than the current average. There are, however, no models which combine both these features even though this is perfectly acceptable in economic situations.
We proposed a new two tier model of learning. In this model, active or ontogenetic’ learning by individuals takes place on a fast time scale and passive or phylogenetic’ learning by society as a whole on a slow time scale, each affecting the evolution of the other. The former is modelled by the Monte Carlo dynamics of statistical physics where the individual updates his strategy to pick moves that yield higher payoffs with higher probability, retaining at the same time `expensive’ moves, albeit with a low probability. The latter may be attributed to random exploration or errors. We superimpose on this, on a slower time scale, the replicator dynamics, describing the evolution of a conglomeration of economic agents under selection pressure from the environment (which may include other such conglomerations). The two dynamics interact through the payoff matrices that enter their velocity fields. This leads to a two time scale evolution, which can then be cast in the framework of singularly perturbed ordinary differential equations. The Monte Carlo dynamics for ontogenetic learning was analyzed in some detail, particularly for the n=2 symmetric case. We carried out numerical simulations which showed a rich variety of qualitative behaviour. In some cases, the two dynamics counter act each other so as to have the system poised in between two distinct kinds of equilibria, reminiscent of (but distinct from) self-organized criticality. In other cases, under similar situations, one obtains sustained oscillations between two equilibria, a phenomenon reminiscent of (but distinct from) noise-induced transitions on one hand and heteroclinic cycles on the other.
We studied the behaviour of the moments of a particle distribution as it is transported through a Hamiltonian system. Functions of moments that remain invariant for an arbitrary linear Hamiltonian system were constructed using Lie algebraic techniques. Several new invariants were obtained. We have also constructed dynamic moment invariants for nonlinear Hamiltonian systems.
We analyzed multi-proxy climatic records from the 74KL marine core and compare it with the oxygen isotope record from the Guliya ice core. These records have important implications for the evolution of the Southwest Monsoon System. We observed three distinct climatic events (viz. 19 ka to 13.5 ka, 13.5 ka to 10 ka and 3.4 ka to present) from the last glacial phase to the present. Even though the first two events are well documented in the literature, our combined analysis of data from both the cores suggested alternative mechanisms which could have complemented and strengthened these events. The third, most recent event (occurring between 3.4 ka to present) does not seem to have been well documented earlier. We proposed a possible mechanism to explain this event. We identified representative elements which better capture the above events as compared to more commonly used elements.
We studied paleoclimatic records from various sites spread around the earth, focusing on the start of the last glacial-interglacial transition. The warming, as recorded in the oxygen isotope record started first in the tropics, then propagated to the Antarctic and then finally to the Arctic. Our analysis of the data suggested that it took about 7.6 ka for onset of climate change to propagate globally. We proposed that the tropical Pacific played a major role in initiating the warming in the tropics. We discussed mechanisms that could have transported this heat from the tropics to Antarctica and then to the Arctic during transition to the interglacial.
We proposed a novel type of free-electron laser (called the standing-wave free-electron laser) for use as a power source in generating high accelerating gradients in next generation accelerators. Theoretical behaviour of this laser was studied. Stability of the system when perturbed was investigated both analytically and numerically. These studies established that the proposed free-electron laser was stable and could be built as a practical power source.
We proposed a novel type of FEL as a practical device to generate brief, intense pulses of radiation. This FEL (which we called the multi-cavity FEL) has many optical cavities where the cavities communicate with each other, not optically, but through the electron beam. The idea is to simply make the FEL optical cavities sufficiently short that the slippage length, in any one optical cavity, is less than the electron pulse width. When electrons reach the end of one cavity they go into the next, but the radiation remains trapped within that cavity. We also studied a mathematical model of this FEL both analytically and numerically. In the one-dimensional approximation, a general expression for the growth rate in the exponential (high current) regime was obtained. In the regime where lethargy is important, expressions were given in the two opposite limits of small and large number of cavities and bunches. Extensive numerical investigations of the nonlinear effects were performed. The multi-cavity FEL was shown to be useful to avoid slippage phenomena, and thus to make pico-second pulses of infra-red radiation. Some examples of this application were studied.
In the wigglers of future FELs, the electron beam will be required to travel over a length of 10m or more in pipes with small diameters. Transverse resistive wall effects could lead to beam breakup during this transport. To investigate this possibility, we solved the equation of motion for a bunched beam analytically. We showed that a steady state solution is reached for times larger than the diffusion time. This solution can either oscillate or grow exponentially with the length of the pipe, depending on the relative magnitudes of the resistive wall effect and the focusing force in the wiggler. The possibility of a significant transient was also studied in great detail.