CAMS Technical Report Series 2005-2006

REPORT 0506-1: Proprioceptor Regulation of Motor Circuit Activity by Presynaptic Inhibition of a Modulatory Projection Neuron

Mark P. Beenhakker¹, Nicholas D. DeLong¹, Shari R. Saideman¹, Farzan Nadim², and Michael P. Nusbaum¹

¹Department of Neuroscience, University of Pennsylvania School of Medicine, 215 Stemmler Hall, Philadelphia, PA 19104

²Department of Mathematical Sciences, New Jersey Institute of Technology & Department of Biological Sciences, Rutgers University, Newark, NJ 07102

Phasically active sensory systems commonly influence rhythmic motor activity via synaptic actions on the relevant circuit and/or motor neurons. Using the crab stomatogastric nervous system (STNS), we have identified a distinct synaptic action by which an identified proprioceptor, the gastro/pyloric muscle stretch receptor (GPR) neuron, regulates the gastric mill (chewing) motor rhythm. Previous work showed that rhythmically stimulating GPR in a gastric mill-like pattern, in the isolated STNS, elicits the gastric mill rhythm via its activation of two identified projection neurons, MCN1 and CPN2, in the commissural ganglia. Here, we determine how activation of GPR with a behaviorally appropriate pattern (active during each gastric mill retractor phase) influences an ongoing gastric mill rhythm via actions in the stomatogastric ganglion, where the gastric mill circuit is located. Stimulating GPR during each retractor phase selectively prolongs that phase and thereby slows the ongoing rhythm. This selective action on the retractor phase results from two distinct GPR actions. First, GPR presynaptically inhibits the axon terminals of MCN1, reducing MCN1 excitation of all gastric mill neurons. Second, GPR directly excites the retractor phase neurons. Because MCN1 transmitter release occurs during each retractor phase, these parallel GPR actions selectively reduce the buildup of excitatory drive to the protractor phase neurons, delaying each protractor burst. Thus, rhythmic proprioceptor feedback to a motor circuit can result from a global reduction in excitatory drive to that circuit, via presynaptic inhibition, coupled with a phase-specific excitatory input that prolongs the excited phase by delaying onset of the subsequent phase.

REPORT 0506-2: Intrinsic Properties, Not Synaptic Dynamics, Determine Proper Phase of Activity in a Central Pattern Generator

Pascale Rabbah¹and Farzan Nadim²

²Department of Mathematical Sciences, New Jersey Institute of Technology and Department of Biological Sciences, Rutgers University, Newark, NJ 07102

Neurons that are involved in the generation of rhythmic motor activity often maintain a relatively constant activity phase despite changes in the rhythm frequency. This maintenance of relative phase is thought to arise from the interplay between the intrinsic properties and the temporal synaptic dynamics among the network neurons. In the rhythmically active pyloric network of the lobster Panulirus interruptus, synaptic connections from the pacemaker ensemble to the follower neurons (lateral pyloric LP and pyloric constrictor PY) are thought to be largely responsible for the proper phase of activity (pacemaker-LP-PY) across all frequencies (0.5-2 Hz) of the pyloric rhythm. We test this hypothesis by characterizing the synapses from the pacemaker ensemble to the LP and PY neurons. Paired comparisons show that these two synapses are not significantly different in strength or in the extent of short-term depression. To examine the level in which intrinsic properties of the follower neurons determine their relative activity phase, we block all chemical synapses within the network and drive the LP and PY neurons rhythmically using artificial synaptic currents with identical strength and dynamics implemented with the dynamic clamp technique. In response to these identical synaptic inputs, the LP and PY neurons maintain the proper relative phase of activity. These results strongly indicate that the relative phase of activity among these follower neurons within the pyloric network is not dictated by their synaptic inputs but is solely determined by their distinct intrinsic properties.

REPORT 0506-3: Convergence of Neuronal Activity Phenotype Result from Spontaneous or Induced Activity via a Common Ionic Mechanism in Adult Isolated Neurons

Rodolfo Haedo¹ and Jorge Golowasch²

¹Biological Sciences, Rutgers University, Newark, NJ 07102, 973-353-5080, roculfo@aol.com

²Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102, 973-353-1267, jorge.p.golowasch@njit.edu

Neurons exhibit long-term changes in excitability in response to a variety of inputs and perturbations. This plasticity of intrinsic properties is necessary for maintaining proper cell and network activity. The adult crustacean pyloric neuronal network can slowly recover rhythmic activity after complete shutdown resulting from permanent removal of neuromodulatory inputs. We use dissociated stomatogastric ganglion (STG) neurons of the crab Cancer borealis as models to study the mechanisms underlying this process. As observed in a different species, STG neurons spontaneously develop a preferred oscillatory activity pattern via gradual changes in excitability, and rhythmic electrical stimulation can regress oscillatory patterns to less excitable states in some cells. However, we show that rhythmic stimulation can more commonly accelerate the emergence of stable oscillatory patterns in cultured crab STG neurons. We find that both spontaneous and activity-induced oscillations correlate with modifications of the same two ionic currents: a Ca++ current increase and a high-threshold K+ current decrease. Dynamic-clamp experiments confirm that these conductance modifications can explain the observed activity-induced changes. We conclude that the majority of stomatogastric ganglion neurons can become endogenous oscillators and retain this capability into adulthood. However, synaptic interactions or trophic factors may prevent the expression of these oscillations in the intact network. The spontaneous excitability changes observed in isolated neurons may explain the spontaneous recovery of rhythmic activity in the functional network after permanent removal of neuromodulatory input in situ.

¹Courant Institute of Mathematical Sciences, New York University, 251 Mercer Street, New York, NY 10012

²Department of Mathematical Sciences, NJIT, 323 Martin Luther King, Jr., Blvd, Newark, NJ 07102

We present a detailed theoretical framework for statistical descriptions of neuronal networks and derive (1+1)-dimensional kinetic equations, without introducing any new parameters, directly from conductance-based integrate-and-fire neuronal networks. We describe the details of derivation of our kinetic equation, proceeding from the simplest case of one excitatory neuron, to coupled networks of purely excitatory neurons, to coupled networks consisting of both excitatory and inhibitory neurons. The dimension reduction in our theory is achieved via novel moment closures. We also describe the limiting forms of our kinetic theory in various limits, such as the limit of mean-driven dynamics and the limit of infinitely fast conductances. We establish accuracy of our kinetic theory by comparing its prediction with the full simulations of the original point-neuron networks. We emphasize that our kinetic theory is dynamically accurate, i.e., it captures very well the instantaneous statistical properties of neuronal networks under time-inhomogeneous inputs.

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

²Department of Mathematics and Programs in Neuroscience and Molecular Biophysics, Florida State University, Tallahassee, FL 32306

³Laboratory for Biological Modeling, NIDDK, National Institutes of Health, Bethesda, MD 20892

Synaptic facilitation is a transient stimulation-induced increase in synaptic transmission strength, a ubiquitous form of short-term synaptic plasticity that may play a role in regulating the activity of synaptically coupled neuronal populations on fast time scales. In their pioneering work, Katz and Miledi (1968) and Rahamimoff (1968) demonstrated the dependence of facilitation on presynaptic calcium influx, and proposed that facilitation results from the accumulation of residual calcium bound to vesicle release triggers. However, in its pure form this bound calcium hypothesis contradicts the evidence that facilitation is reduced by exogenous calcium buffers, which suggests the importance of free calcium for this form of synaptic plasticity. This led to a widely-held view that facilitation must depend solely on the accumulation of calcium in free form. Here we consider a more realistic implementation of the bound calcium mechanism, taking into account spatial diffusion of calcium ions, and show that a model with a slow calcium unbinding step can retain sensitivity to free residual calcium.

We demonstrate that this hybrid free/bound calcium model agrees with the facilitation accumulation time course exhibited by the crayfish inhibitor neuromuscular junction (NMJ), and relies on fewer assumptions than the most recent variations of the free residual calcium hypothesis.

Further, we show that the hybrid mechanism is consistent with the experimental results of Kamiya and Zucker (1994), which revealed that a photolytic liberation of a fast calcium buffer decreases the synaptic response within milliseconds. We conclude that a calcium binding process with a slow unbinding time step (10s of ms) is a viable mechanism of synaptic facilitation at some synapses, and discuss the experimental evidence for such a mechanism.

REPORT 0506-6: Quenching and Propagation of Combustion Fronts in Porous Media

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

In this short note we study the model of subsonic detonation introduced by Sivashinsky. The model is described by the system of reaction-diffusion equations involving temperature, pressure and concentration of deficient reactant. It is shown that initial data with small support lead to the quenching (decay of solution). In contrast initial data with support large enough lead to propagation with finite velocity.

REPORT 0506-7: From Immersed Boundary Method to Immersed Continuum Method

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The objective of this paper is to present an overview of the newly proposed immersed continuum method in conjunction with the traditional treatment of fluid-structure interaction problems, the immersed boundary method, the extended immersed boundary method, the immersed finite element method, and the fictitious domain method. In particular, the key aspects of the immersed continuum method in comparison with the immersed boundary method are discussed. The immersed continuum method retains the same strategies employed in the extended immersed boundary method and the immersed finite element method, namely, the independent solid mesh moves on top of a fixed or prescribed background fluid mesh, and employs fully implicit time integration with a matrix-free combination of Newton-Raphson and GMRES iterative solution procedures. Therefore, the immersed con-tinuum method is capable of handling compressible fluid interacting with compressible solid. Several numerical examples are also presented to demonstrate that the proposed immersed continuum method is a good candidate for multi-scale and multi-physics modeling platform.

REPORT 0506-8: Immersed Finite Element Method and Its Applications to Biological Systems

Wing Kam Liu¹, Yaling Liu¹, David Farrell¹, Lucy Zhang¹, X. Sheldon Wang², and Others

¹Department of Mechanical Engineering, Northwestern University, Evanston, IL 60208

²Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

REPORT 0506-9: Using Experimental Methods and Mathematical Modeling to Determine Gap Junction Coupling in Neural Networks

Diana Martinez, Matt Malej, Angelie Mascarinas, Farzan Nadim, and Jorge Golowasch

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The existence of electrical synapses between neurons has been determined to play an important role in communication and network synchronization. The location of the gap junctional communication between neurons is of functional importance because of the complex neuronal morphology and the non-uniform distribution of voltage-gated ionic channels. However, accurately locating the gap junctions has been an elusive task. We performed intracellular recordings of two gap-junctionally coupled neurons in the crab Cancer borealis and examined signal transfer between these neurons using a variety of stimuli. We constructed a simplified mathematical model of the neurons to predict the location of the gap junctional coupling using the obtained experimental data.

REPORT 0506-10: Border Properties of Allopatric and Sympatric Plant Species Interactions

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

We investigated the borders between 15 different combinations of perennial plant species and determined that sympatric species pairs (species that originate from similar geographic areas) have borders with less overlap between them than those between allopatric species (species that originate from different areas). We used a grid transect method to assess/ /aboveground interactions between 9 pairs of sympatric species and 6 pairs of allopatric species. Our experimental results were then matched to a general mathematical model, which reproduces the observations witnessed in the field and aims at separating competitive exclusion and active interference. It is discussed how to utilize the model to make qualitative and quantitative predictions about the ecology between two competing species.

²Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

We address the question of whether or not a periodic train of excitatory synaptic inputs recruits an excitable cell, such that it fires repeatedly, or does not recruit a cell, such that it fails to fire, possibly after some transient. In particular, we study the scenarios of one or two inputs per period; in the latter case, the degree of synchrony of the inputs is a crucial factor in recruitment. We establish rigorous geometric conditions that pinpoint the transition between recruitment and non-recruitment as the degree of synchrony between input pairs, or other input parameters, is varied. These conditions can be used to determine whether a particular temporal relation between periodic input pairs leads to recruitment or not and to prove, in certain parameter regimes, that recruitment can only occur when the inputs are sufficiently closely synchronized. The concepts in this paper are derived for both the integrate-and-fire neuron and the theta neuron models. In the former, the location in phase space of the unique fixed point of a relevant two-dimensional map determines firing, while in the latter, it is the existence or lack of existence of a fixed point of the map which does so. These results are discussed in the context of recruitment of cells into localized activity patterns.

REPORT 0506-12: Properties of Fast and Slow Endogenous Buffers at the Presynaptic Terminals of the Crayfish Neuromuscular Junction

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, New Jersey 07102

Endogenous Ca2+ buffers (ECBs) typically bind 98% to 99.9% of cytosolic Ca2+, and therefore play an important role in regulating synaptic transmission and other Ca2+-dependent processes; however, the biophysical characteristics of buffers are unknown in most cell types.

We attempt to infer the kinetic properties of ECBs at the crayfish inhibitor neuromuscular junction, exploring the recent suggestion by Lin et al. (2005) that an endogenous buffering system consisting of one fast and one slow buffer can account for the time course of Ca2+ transients recorded in this preparation. It was proposed that the ECB with slow Ca2+ binding-unbinding kinetics controls the decay time course of the AP-evoked Ca2+ transient (~ 50 ms), whereas the fast buffer reduces the amplitude of the transient, and accounts for the low degree of saturation of the indicator dye during the recording. To estimate the properties of these two buffers, we compare the time course of fluorescence transient evoked by two broadened action potentials to the one obtained through computational modeling, both in the control condition and upon application of EGTA. Our results indicate that the fast buffer comprises about 10 to 15% of the total buffering capacity, previously estimated by Tank et al. (1995) to equal about 600 in this preparation, with the remaining 85-90% of the buffering ratio contributed by the slow buffer. We infer that the slow buffer has low affinity (>1 uM), is present in large concentration (>500 uM), and we estimate its unbinding rate to range between 1/sec and 4/sec. Finally, we also determine the bounds on the affinity (>0.4 uM) and the total concentration (>20 uM) of the fast ECB.

REPORT 0506-13: Generalized Poincare-Bertrand Formula on a Regular Surface in Three Dimensions

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

²Department of Mathematics, Princeton University, Fine Hall, Washington Road Princeton, NJ 08544

The Poincare-Bertrand formula is generalized to the case of singular integrals on a regular surface in three dimensions. The generalized formula is a natural extension of the original Poincare-Bertrand formula concerning about two repeated Cauchy's principal integrals on a sufficiently smooth curve in two dimensions. The formula is expected to be useful in constructing a second kind integral equation formulation for scattering by open surfaces in three dimensions in the acoustic and electromagnetic environments.

REPORT 0506-14: Scaling Law for Second-Order Hyperpolarizability of Periodic Finite Chain Under Su-Schrieffer-Heeger Model

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The analytical format of the second hyperpolarizability with N double-bonds finite chain is obtained under the Su-Schrieffer-Heeger model of trans-polyacetylene. The feature of the second hyperpolatizability versus N can be qualitatively compared with the existing experiments. Our results provide the physical understanding of the main feature of saturation behaviors in polyenes. It suggests that more physical effects should be included to do an accurate quantitative comparison with experiments.

REPORT 0506-15: Surface Tension in Incompressible Rayleigh-Taylor Mixing Flow

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

We study the effect of surface tension on the incompressible Rayleigh-Taylor instability. We modify Goncharov's local analysis [1] to consider the surface tension effect on the Rayleigh-Taylor bubble velocity. The surface tension damps the linear instability and reduces the nonlinear terminal bubble velocity. We summarize the development of a finite-volume, particle-level-set, two-phase flow solver with an adaptive Cartesian mesh, and results from convergence and validation studies of this two-phase flow solver are provided. We use this code to simulate the single-mode, viscous Rayleigh-Taylor instability with surface tension, and good agreement in terminal bubble velocity is found when compared with analytic results. We also simulate the immiscible Rayleigh-Taylor instability with random initial perturbations. The ensuing mixing flow is characterized by the effective mixing rate and the flow anisotropy. Surface tension tends to reduce the effective mixing rate and homogenizes the Rayleigh-Taylor mixing flow. Finally we provide a scaling argument for detecting the onset of the quadratic, self-similar Rayleigh-Taylor growth.

REPORT 0506-16: Limits of the Potential Flow Approach to the Single-Mode Rayleigh-Taylor Problem

P. Ramaprabhu¹, Guy Dimonte¹, Yuan-Nan Young², Alan C. Calder^3, and Bruce Fryxell⁴

¹Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545

²Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

⁴Aerospace Engineering Department, Georgia Institute of Technology, Atlanta, GA 30332

We report on the behavior of a single-wavelength Rayleigh-Taylor flow at late-times. The calculations were performed in a long square duct (λxλx8λ), using four different numerical simulations. In contradicitn with potential flow theories that predict a constant terminal velocity, the single-wavelength Rayleigh-Taylor problem exhibits late-time acceleration. The onset of acceleration occurs as the bubble penetration depth exceeds the diameter of bubbles, and is observed for low and moderate density differences. We provide a phenomenological description of the observed acceleration, and ascribe this behavior to the formation of Kelvin-Helmholtz vortices on the bubble-spike interface that diminish the friction drag. For large density ratios, the formation of secondary instabilities is suppressed, and the bubbles remain terminal consistent with potential flow models.

Glass microelectrodes are used widely in experimental studies of the electrophysiology of biological cells and their membranes. However, the pulling of these electrodes remains an art, based on trial-and-error. Following Huang et al. (2003), we derive a one-dimensional model for the stretching of a hollow glass tube that is being radiatively heated. Our framework allows us to consider two commonly used puller designs, that is, horizontal (constant force) and vertical (variable force) pullers. We derive explicit solutions and use these solutions to identify the principal factors that control the final shape of the microelectrodes. The design implications for pullers also are discussed.

REPORT 0506-18: Traveling Waves in Coupled Reaction-Diffusion Models with Degenerate Sources

We consider a general system of coupled nonlinear diffusion equations that are characterized by having degenerate source terms and thereby not having isolated rest states. Using a general form of physically relevant source terms, we derive conditions that are required to trigger traveling waves when a stable uniform steady-state solution is perturbed by a highly localized disturbance. We show that the degeneracy in the source terms implies that traveling waves have a number of surprising properties that are not present for systems with non-degenerate source terms. We also show that such systems can lead to a pair of waves that initially propagate outwards from the disturbance, slow down, and reverse direction before ultimately colliding and annihilating each other.

REPORT 0506-19: Log-linear Modeling Under Generalized Inverse Sampling Scheme

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

This paper discusses the log-linear model for multi-way contingency table, where the cell values represent the frequency counts that follow an extended negative multinomial distribution. This is an extension of negative multinomial log-linear model described by Evans (1989). The parameters of the new model are estimated by maximum likelihood method. The likelihood ratio test for the general log-linear hypothesis is also derived. A practical application of the log-linear model under the generalized inverse sampling scheme has also been demonstrated by an example.

REPORT 0506-20: Steady Deformation and Tip-Streaming of a Slender Bubble with Surfactant in an Extensional Flow

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

Slender body theory is used to investigate the steady-state deformation and time-dependent evolution of an inviscid, axisymmetric bubble in zero-Reynolds-number extensional flow, when insoluble surfactant is present on the bubble surface. The asymptotic solutions reveal steady ellipsoidal bubbles covered with surfactant, and, at increasing deformation, solutions distinguished by a cylindrical, surfactant-free central part, with stagnant surfactant caps at the bubble end-points. The bubble shapes are rounded near the end-points, in contrast to the pointed shapes (Buckmaster 1972) found for clean, inviscid bubbles. Simple expressions are derived relating the capillary number Q to the steady bubble slenderness ratio. These show that there is a critical value of capillary number Qc above which steady solutions no longer exist. Equations governing the time-evolution of a slender inviscid bubble with surfactant, valid for large capillary number, are also derived. Numerical solutions of the slender bubble equations for Q greater than Qc exhibit spindle shapes with tip-streaming filaments.

REPORT 0506-21: Mathematical Model for the Rupture of Cerebral Saccular Aneurysms through Three-dimensional Stress Distribution in the Aneurysm Wall

Hans R Chaudhry^1, Dawn A. Lott², Charles J. Prestigiacomo³, and Thomas W. Findley⁴

¹Department of Biomedical Engineering, New Jersey Institute of Technology, Newark, New Jersey and War-Related Illness and Injury Study Center, VA Medical Center, East Orange, New Jersey.

²Departments of Biotechnology and Mathematics, Applied Mathematics Research Center, Delaware State University, Dover, Delaware and Center for Applied Mathematics and Statistics, New Jersey Institute of Technology, Newark, New Jersey

³Department of Neurological Surgery, New Jersey Medical School, University of Medicine and Dentistry of New Jersey, Newark, New Jersey and Department of Biomedical Engineering, New Jersey Institute of Technology, Newark, New Jersey

⁴War-Related Illness and Injury Study Center, VA Medical Center, East Orange, New Jersey

A mathematical model for the rupture of cerebral saccular aneurysms is developed through the analysis of three-dimensional stress distribution in the aneurysm wall. We assume in this paper, that a saccular aneurysm resembles a thin spherical shell (a spherical membrane), and then develop a strain energy function valid for finite strain to analyze 3-dimensional stress distribution in the aneurysm wall. We find that rupture occurs when the ratio of the wall thickness to the radius of the aneurysm is 6.1 X 10-3. We also conclude from our analysis that rupture can occur when the ratio of thickness to radius of the parent aneurysm equals the ratio of thickness to radius of the daughter aneurysm. These findings may be useful to the neurosurgeon to help predict the rupture potential in patients presenting with unruptured aneurysms.

REPORT 0506-22: Maximum a Posteriori Multiple Source Localization with Gibbs Sampling

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

Multiple source localization in underwater environments is approached within a matched-field processing framework. A Maximum a Posteriori Estimation method is proposed that estimates source location and spectral characteristics of multiple sources via Gibbs Sampling. The method facilitates localization of weak sources which are typically masked by the presence of strong interferers. A performance evaluation study based on Monte Carlo simulations shows that the proposed Maximum a Posteriori Estimation approach is superior to simple coherent matched-field interference cancellation. The proposed method is also tested on the estimation of the number of sources present, providing probability distributions in addition to point estimates for the number of sources.

REPORT 0506-23: An explicit formulation for the evolution of nonlinear surface waves interacting with a submerged body

²Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

An explicit formulation to study nonlinear waves interacting with a submerged body in an ideal fluid of infinite depth is presented. The formulation allows one to decompose the nonlinear wave-body interaction problem into the body and free-surface problems. After the decomposition, the body problem satisifies a modified body boundary condition in an unbounded °uid domain, while the free-surface problem satistisfies modified nonlinear free-surface boundary conditions. It is then shown that the nonlinear free-surface problem can be further reduced to a closed system of two nonlinear evolution equations expanded in infinite series for the free-surface elevation and the velocity potential at the free surface. For numerical experiments, the body problem is solved using a distribution of singularities along the body surface and the system of evolution equations truncated at third order in wave steepness is then solved using a pseudo-spectral method based on the Fast Fourier Transform. A circular cylinder translating steadily near the free surface is considered and it is found that our numerical solutions show excellent agreement with the fully nonlinear solution using a boundary integral method. We further validate our solutions for a submerged circular cylinder oscillating vertically or fixed under incoming nonlinear waves with other analytical and numerical results.

REPORT 0506-24: The Effect of a Background Shear Current on Large Amplitude Internal Solitary Waves

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

Large amplitude internal waves interacting with a linear shear current in a system of two layers of different densities are studied using a set of nonlinear evolution equations derived under the long wave approximation without the smallness assumption on the wave amplitude. For the case of uniform vorticity, solitary wave solutions are obtained under the Boussinesq assumption for a small density jump, and the explicit relationship between the wave speed and the wave amplitude is found. It is shown that a linear shear current modifies not only the wave speed, but also the wave profile drastically. For the case of negative vorticity, when compared with the irrotational case, a solitary wave of depression traveling in the positive x-direction is found to be smaller, wider, and slower, while the opposite is true when traveling in the negative x-direction. In particular, when the amplitude of the solitary wave propagating in the negative x-direction is greater than the critical value, a stationary recirculating eddy appears at the wave crest.

REPORT 0506-25: Breeding Birds on Small Islands: Island Biogeography or Optimal Foraging?

Gareth J. Russell¹, Jared M. Diamond², Timothy M. Reed³, and Stuart L. Pimm⁴

¹Department of Mathematical Sciences, Indiana University Purdue University, Indianapolis, IN 46202, USA

²Department of Mathematical Sciences, New Jersey Institute of Technology, University Heights, Newark, NJ 07102 USA

³Department of Mathematical Sciences, Tel Aviv University, Tel Aviv 69978, Israel

The classical relation between the flame speed and the stretch, employed in modeling flame-flow interaction, is valid only for positive Markstein lengths (high Lewis numbers). At negative Markstein lengths (low Lewis numbers) the corresponding dynamical system suffers short-wavelength instability, making the associated initial value problem ill-posed. In this study the difficulty is resolved by incorporation of higher-order effects. As a result one ends up with a reduced model based on a coupled system of second-order dynamic equations for the flame interface and its temperature. As an illustration the new model is applied for description of diffusively unstable stagnation-point flow flames.

REPORT 0506-27: A Mathematical Modelling Approach to One-Day Cricket Batting Orders

¹School of Mathematical Sciences, Faculty of Science, Monash University, Australia

²Center for Applied Mathematics and Statistics, Department of Mathematical Sciences, New Jersey Institute of Technology Newark, NJ, USA 07102

While scoring strategies and player performance in cricket have been studied, there has been little published work about the influence of batting order with respect to One-Day cricket. We apply a mathematical modelling approach to compute efficiently the expected performance (runs distribution) of a cricket batting order in an innings. Among other applications, our method enables one to solve for the probability of one team beating another or to find the optimal batting order for a set of 11 players. The influence of defence and bowling ability can be taken into account in a straightforward manner.

In this presentation, we outline how we develop our Markov Chain approach to studying the progress of runs for a batting order of non-identical players along the lines of work in baseball modelling by Bukiet et al. (1997). We describe the issues that arise in applying such methods to cricket, discuss ideas for addressing these difficulties and note limitations on modelling batting order for One-Day cricket.

By performing our analysis on a selected subset of the possible batting orders, we apply the model to quantify the influence of batting order in a game of One Day cricket using available real-world data for current players.

REPORT 0506-28: Kink-antikink collisions in the phi-four equation: The n-bounce resonance and the separatrix map

¹Department of Mathematical Sciences, New Jersey Institute of Technology, University Heights, Newark, NJ 07102

We provide a detailed mathematical explanation of a phenomenon known as the two-bounce resonance observed in collisions between kink and anti-kink traveling waves of the phi-four equations of mathematical physics. This behavior was discovered numerically in the 1980's by Campbell and his collaborators and subsequently discovered in several other equations supporting traveling waves. We first demonstrate the effect with new high-resolution numerical simulations. A pair of kink- like traveling waves may coalesce into a localized bound state or may reflect off each other. In the two bounce-resonance, they first coalesce, but later escape each others' embrace, with a very regular pattern governing the behaviors. Studying a finite-dimensional ``collective coordinates'' model, we use geometric phase-plane based reasoning and matched asymptotics to explain the mechanism underlying the phenomenon, including the origin of several mathematical assumptions needed by previous researchers. We derive a separatrix map for this problem---a simple algebraic recursion formula that explains the complex fractal-like dependence on initial velocity for kink-antikink interactions.

REPORT 0506-29: Long-Wave Linear Stability Theory for Two-Fluid Channel Flow Including Compressibility Effects

¹Department of Chemical and Materials Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G6

²Department of Mathematical Sciences, Center for Applied Mathematics and Statistics, New Jersey Institute of Technology Newark, NJ 07102

We present the linear stability of the laminar flow of an immiscible system of a compressible gas and incompressible liquid separated by an interface with large surface tension in a thin inclined channel. The flow is driven by an applied pressure drop and gravity. Following the air-water case, which is found in a variety of engineering systems, the ratio of the characteristic values of the gas and liquid densities and viscosities are assumed to be disparate. Under lubrication approximation, and assuming ideal gas behavior and isothermal conditions, this approach leads to a coupled nonlinear system of partial differential equations describing the evolution of the interface between the gas and the liquid and the streamwise density distribution of the gas. This system also includes the effects of viscosity stratification, inertia, shear, and capillarity. A linear stability analysis that allows for physically relevant nonzero pressure-drop base state is then performed. In contrast to zero-pressure drop case which is amenable to the classical normal-mode approach, this configuration requires solving numerically a boundary-value problem for the gas density and interfacial deviations from the base state in the streamwise coordinate.

We find that the effect of gas compressibility on the interfacial stability in the limit of vanishingly small wavenumber is destabilizing, even for Stokes flow in the liquid. However, for finite wavenumber disturbances, compressibility may have stabilizing effects. In this regime, sufficient shear is required to destabilize the flow.

REPORT 0506-30: Stabilization of Nonlinear Velocity Profiles in Athermal Systems Undergoing Planar Shear Flow

¹Department of Mechanical Engineering, Yale University, New Haven, CT 06520-8284

²Department of Mechanical Engineering and Department of Physics, Yale University, New Haven, CT 06520-8120

³Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

We perform molecular dynamics simulations of model granular systems undergoing boundary-driven planar shear flow in two spatial dimensions to develop a more complete understanding of how dense particulate systems respond to applied shear. In particular, we are interested in determining when these systems will possess linear velocity profiles and when they will develop highly localized velocity profiles in response to shear. We showed in previous work on similar systems that nonlinear velocity profiles form when the speed of the shearing boundary exceeds the speed of shear waves in the material. However, we find that nonlinear velocity profiles in these systems are unstable at very long times. The degree of nonlinearity slowly decreases in time; the velocity profiles become linear when the granular temperature and density profiles are uniform across the system at long times. We measure the time t_l required for the velocity profiles to become linear and find that t_l increases as a power-law with the speed of the shearing boundary and increases rapidly as the packing fraction approaches random close packing. We also performed simulations in which differences in the granular temperature across the system were maintained by vertically vibrating one of the boundaries during shear flow. We find that nonlinear velocity profiles form and are stable at long times if the difference in the granular temperature across the system exceeds a threshold value that is comparable to the glass transition temperature in an equilibrium system at the same average density. Finally, the sheared and vibrated systems form stable shear bands, or highly localized velocity profiles, when the applied shear stress is lowered below the yield stress of the static part of the system.

REPORT 0506-31: Elastic Energy, Fluctuations and Temperature for Granular Materials

¹Department of Mathematical Sciences & Center for Applied Mathematics & Statistics New Jersey Institute of Technology, Newark, NJ 07102

²Department of Physics and Center for Nonlinear and Complex Systems Duke University, Durham NC, 27708-0305

In our recent work (Europhys. Lett. Vol. 67, 205 (2004)) we have shown that in granular systems characterized by large volume fractions, the elastic energy dominates the kinetic energy, and energy fluctuations are primarily elastic in nature. As a logical consequence of this observation, we have started exploring possible generalizations of the concept of granular temperature to dense, jammed systems where kinetic granular temperature is not expected to be relevant, at least from the energetic point of view. Therefore, we have introduced generalized granular temperature, which turns out to be roughly consistent with a temperature based on the equilibrium statistical mechanics. In this paper, we discuss the influence of various system properties on this new generalized granular temperature. These properties include the shearing rate, as well as the material properties such as stiffness, elasticity and friction.

REPORT 0506-32: Numerical Methods for Solving Kinetic Equations of Neuronal Network Dynamics

¹Courant Institute of Mathematical Sciences, New York University, New York, NY 10012

²Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

Recently developed kinetic theory for neuronal network dynamics has been demonstrated to be a powerful theoretical framework for investigating coarse-grained dynamical properties of neuronal networks. The kinetic equations are a system of (1+1)-dimensional nonlinear partial differential equations (PDE) on a bounded domain with the following features: (i) the boundary conditions are nonlinear and they are themselves a functional of the present solution; (ii) the PDEs themselves are self-consistently specified by parameters which are functions of the boundary values of the solution; and (iii) the PDEs can be stiff in space and time. Numerical methods are presented here for efficiently and accurately solving these kinetic equations. The essential ingredients in our numerical methods include (i) the system is discretized in time with an implicit Euler method within a spectral deferred correction framework, therefore, the PDEs of the kinetic theory are reduced to a sequence, in time, of boundary value problems (BVPs) with nonlinear boundary conditions; (ii) a set of auxiliary parameters is introduced to recast the original BVP with nonlinear boundary conditions as BVPs with linear boundary conditions --- with additional algebraic constraints on the auxiliary parameters; and (iii) a careful combination of two Newton's iterate for solving the associated algebraic constraints is constructed to achieve quadratic convergence for obtaining the solutions with self-consistent parameters. It is shown that a simple fixed-point iteration can only achieve a linear convergence for the self-consistent parameters. The practicability and efficiency of our numerical methods for solving the kinetic equations are illustrated with numerical examples. It is further demonstrated that the kinetic equations can very well capture the coarse-grained dynamical properties of integrate-and-fire neuronal networks.

REPORT 0506-33: Soliton Broadening under Random Dispersion Fluctuations: Importance Sampling Based on Low-Dimensional Reductions

¹Department of Mathematical Sciences, New Jersey Institute of Technology, 323 Martin Luther King, Jr. Blvd., Newark, NJ 07102

²Department of Mathematics, The College of Staten Island, 2800 Victory Blvd., Staten Island, NY 10314

³Department of Mathematics, The University of North Carolina at Chapel Hill, CB #3250, Phillips Hall, Chapel Hill, NC 27599

We demonstrate that dispersion-managed solitons are less likely to experience critical broadening under the influence of random dispersion fluctuations than are solitons of the integrable nonlinear Schroedinger equation, and that this robustness increases with map strength from the constant-dispersion(integrable) limit to the large-map-strength limit. To achieve this, we exploit a separation of scales in dispersion-managed soliton dynamics to implement an importance-sampled Monte Carlo approach that determines the probability of rare broadening events. This approach reconstructs the tails (i.e., the regions of practical importance) of probability distribution functions with an efficiency that is several orders of magnitude greater than conventional Monte Carlo simulations. We further show that the variational approach with an appropriately scaled ansatz is surprisingly good at capturing the effect of random dispersion on pulse broadening; where it fails, it can still be used to guide very efficient simulation of the original equation.

REPORT 0506-34: Renormalization Group Reduction of Pulse Dynamics in Thermally Loaded Optical Parametric Oscillators

¹Department of Mathematical Sciences, New Jersey Institute of Technology, 323 Martin Luther King, Jr. Blvd., Newark, NJ 07102

We derive a perturbed parametrically forced nonlinear Schroedinger equation to model pulse evolution in an optical parametric oscillator with absorption-induced heating. We apply both a rigorous renormalization group technique and the formal Wilsonian renormalization group to obtain a low-dimensional system of equations which captures the mutual interaction of pulses as well as their response to the thermally induced potential. We compare the methodologies of the two approaches and find that the two reduced systems agree to leading order, and compare well to simulations of the full equation, predicting the formation of stably bound pulse pairs.

REPORT 0506-35: A Method to Compute Statistics of Large, Noise-Induced Perturbations of Nonlinear Schrödinger Solitons

¹Department of Mathematical Sciences, New Jersey Institute of Technology, 323 Martin Luther King, Jr. Blvd., Newark, NJ 07102

²Department of Mathematics, State University of New York at Buffalo, Buffalo, NY 14260

³Department of Engineering Sciences and Applied Mathematics Northwestern University Evanston, IL 60208

We show in detail how to apply importance sampling to the numerical simulation of large, noise-induced perturbations in soliton-based optical transmission systems governed by the nonlinear Schroedinger equation. The method allows one to concentrate numerical Monte Carlo simulations around the noise realizations that are most likely to produce the large pulse deformations connected with errors, and yields computational speedups of several orders of magnitude over standard Monte Carlo simulations. We demonstrate the method by using it to calculate the probability density functions associated with pulse amplitude, frequency and timing fluctuations in a prototypical soliton-based communication system.

REPORT 0506-36: Neuromodulation of Short-Term Synaptic Dynamics Examined in a Mechanistic Model Based on Kinetics of Calcium Currents

Network plasticity arises in large part due to the effects of exogenous neuromodulators. We investigate the neuromodulatory effects on short-term synaptic dynamics. The synapse from the lateral pyloric (LP) to the pyloric dilator (PD) neuron in the pyloric network of the crab C. borealis has both spike-mediated and non-spike-mediated (graded) components. Previous studies have shown that the graded component of this synapse exhibits short-term depression. Recent results from our lab indicate that in the presence of neuromodulatory peptide proctolin, low-amplitude presynaptic stimuli switch the short-term dynamics of this graded component from depression to facilitation. In this study, we show that this facilitation is correlated with the activation of a presynaptic inward current that is blocked by Mn2+ suggesting that it is a slowly-accumulating Ca2+ current. We modify a mechanistic model of synaptic release by assuming that the low-voltage-activating Ca2+ current in our system is composed of two currents with fast (ICaF) and slow (ICaS) kinetics. We show that if proctolin adjusts the activation rate of ICaS, this leads to accumulation of local intracellular Ca2+ in response to multiple presynaptic voltage stimuli which, in turn, results in synaptic facilitation. Additionally, we assume that proctolin increases the maximal conductances of Ca2+ currents in the model, consistent with the increased synaptic release found in the experiments. We find that these two presynaptic actions of proctolin in the model are sufficient to describe its actions on the short-term dynamics of the LP to PD synapse.

REPORT 0506-37: Combining Synaptic and Cellular Resonance in a Feed-Forward Neuronal Network

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

We derive a mathematical theory to explain the subthreshold resonance response of a neuron to synaptic input. The theory shows how a neuron combines information from its intrinsic resonant properties with those of the synapse to determine the neuron's generalized resonance response. Our results show that the maximal response of a postsynaptic neuron can lie between the preferred intrinsic frequency of the neuron and the synaptic resonance frequency. We compare our theoretical results to parallel findings on experiments of the crab pyloric central pattern generator.

REPORT 0506-38: Modeling Recovery of Rhythmic Activity: Hypothesis for the Role of a Calcium Pump

Yili Zhang¹ and Jorge Golowasch²

¹Federated Department of Biological Sciences, Rutgers University, Newark, NJ

²Department of Mathematical Sciences and Federated Department of Biological Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The pyloric network of crustaceans is a model system for the study of the recovery of function after perturbation/injury of a central pattern generating network. The network is well characterized anatomically and functionally, yet the cellular mechanism underlying the stabi-lization or recovery of its activity is not known. In a previous theoretical study long-term ac-tivity-dependent regulation of ionic conductances was shown to be sufficient to explain the recovery of rhythmic activity after it is lost due to removal of central input. This model, how-ever, did not capture the complex temporal activity dynamics (bouting) that follows decen-tralization and that precedes the final stable recovery. Here we build a model of a conditional pacemaker neuron whose ionic conductance levels depend on activity as before, but also in-cludes a slow activity-dependent regulation of the Ca++ uptake (and release). Intracellular Ca++ sensors, representing enzymatic pathways, regulate the Ca++ pump activity as well as Ca++ and K+ conductances. Our model suggests that the activity-dependent regulation of Ca++ uptake as well as ionic currents interact to generate the complex changes in pyloric activity that follows decentralization. Supported by NIMH 64711 and NSF IBN-0090250.

REPORT 0506-39: Ionic Mechanism Underlying Recovery of Rhythmic Activity in Adult Isolated Neurons

Rodolfo J. Haedo¹ and Jorge Golowasch²

Neurons exhibit long-term excitability changes necessary for maintaining proper cell and network activity in response to various inputs and perturbations. For instance, the adult crustacean pyloric network can spontaneously recover rhythmic activity after complete shutdown resulting from permanent removal of neuromodulatory inputs. Dissociated lobster stomatogastric ganglion (STG) neurons have been shown to also spontaneously recover oscillatory activity via gradual changes in excitability but rhythmic electrical stimulation can eliminate these oscillatory patterns in some cells. The ionic mechanisms underlying these changes are only partially understood. Here we used dissociated crab STG neurons to study the ionic mechanisms underlying the spontaneous recovery of rhythmic activity and stimulation-induced modifications of activity. Similar to lobster neurons, rhythmic activity spontaneously recovers in crab STG neurons. We show that rhythmic stimulation can eliminate, but more commonly accelerate the emergence of, stable oscillatory activity. Our main finding is that up-regulation of a Ca++-current and down-regulation of a high-threshold K+-current is the mechanism underlying both spontaneous and activity-induced oscillations. Moreover, we find no difference in the activity patterns or the underlying ionic current involved between neurons of the fast pyloric and the slow gastric mill networks during the first ten days in isolation. Dynamic-clamp experiments confirm that these conductance modifications can explain the observed activity-induced changes. We conclude that spontaneous and stimulation-induced excitability changes in STG neurons both result in endogenous oscillatory activity via a Ca++-dependent mechanism involving the same two conductances, in both cases irrespective of neuronal identity.

REPORT 0506-40: A Family of Probability Generating Functions Induced by Shock Models

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The literature regarding the non-parametric aging properties are quite big. Bhattacharjee(2005) has introduced a new notion of non-parametric aging property known as Strong decreasing Failure rate (SDFR). This paper explores the necessary and sufficient condition for which this nonparametric aging property is preserved under EMP shock model. It has been proved that the discrete SDFR property is transmitted to continuous version of SDFR under shock model operation. An example is illustrated to show that the converse is not true.

REPORT 0506-41: Analysis of Biological Neurons Via Modeling And Rule Mining

Tomasz G. Smolinski¹, Pascale Rabbah^2,3, Cristina Soto-Treviño², Farzan Nadim ^2,3, and Astrid A. Prinz¹

Due to experimental constraints, measurement errors and variability, analyzing how the activity of biological neurons depends on cellular parameters can be difficult. Computational modeling of neurons allows for exploration of many parameter combinations and various types of neuronal activity, without requiring a prohibitively large number of "wet" experiments. Databases of model neurons created through parameter exploration can, however, be very extensive. There thus is a need for an automated analysis of high-dimensional parameter spaces to explain how neuronal parameters influence the output activity of the modeled cells. In this article, we propose an evolutionary algorithms-based pseudoassociation rule mining methodology to deal with this task.

REPORT 0506-42: A Global Description of Solutions to Nonlinear Perturbations of the Wiener-Hopf Integral Equations

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

²Department Mathematics and Statistics, York University, Toronto, Ontario, Canada M3J 1P3

³Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102, USA

REPORT 0506-44: Spatial Buffering Mechanism: Mathematical Model and Computer Simulations

¹Institute of Medical Science, University of Toronto, Toronto, Ontario, M5S 1A8, Canada

²Pacific Institute for the Mathematical Sciences, Vancouver, BC, V6T 1Z2 Canada

³Department of Mathematics and Statistics, York University, 4700 Keele Street, Toronto, ON M3J 1P3, Canada

⁴Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102, USA

It is generally accepted that the spatial buffering mechanism is important to buffer extracellular-space potassium in the brain-cell microenvironment. In the past, this phenomenon, generally associated with glial cells, has been treated analytically and numerically using a simplified one-dimensional description. The present study extends the previous research by using a novel numerical scheme for the analysis of potassium buffering mechanisms in the extracellular brain-cell microenvironment. In particular, a lattice-cellular automaton was employed to simulate a detailed two-compartment model of a two-dimensional brain-cell system. With this numerical approach, the present study elaborates upon previous theoretical work on spatial buffering (SB) by incorporating a more realistic structure of the brain cell microenvironment, which was not feasible earlier. We use the experimental paradigm consisting of iontophoretic injection of KCl to study the SB mechanism. Our simulations confirmed the results reported in the literature obtained by an averaged model. The results also show that the additional effects captured by a simplified twodimensional geometry do not alter significantly the conclusions obtained from the averaged model. The details of applying such a numerical method to the study of ion movements in cellular environments, as well as its potential for future study, are discussed.

¹Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

A model of experiments where electrical stimulation of the eye of human subjects results in the perception of evenly spaced lines, or phosphenes, is presented. The model is a two-dimensional grid of integrate and fire oscillators that captures the important experimental characteristics of line-creation when a sinusoidal current injection is used. The spatio-temporal behavior of the lines, once formed, is also reproduced. A reduced model consisting of an evolution/convolution equation on the real line is analyzed, and it is shown that stationary solutions with arbitrarily located discontinuities exist and are linearly stable. Traveling waves are numerically shown to exist when the coupling is both sufficiently strong and biased which accounts for the movement of the lines in the experiments.

REPORT 0506-46: On the Calculation of Convolutions with Gaussian Kernels

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The calculation of convolutions with Gaussian kernels arises in many applications. The accuracy of a family of recently developed numerical integration rules is studied in this paper. In spite of the attractive implementation properties of the method, the poor convergence properties greatly restrict the situations in which the method should be used.

REPORT 0506-47: Mesoscopic Simulation for Self-Organization in Surface Processes

Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

The self-organization of particles in a system through a diffusive mechanism is known as Ostwald ripening. This phenomenon is an example of a multiscale problem in that the microscopic level interaction of the particles can greatly impact the macroscale or observable morphology of the system. The mesoscopic model of this physical situation is a stochastic partial differential equation which can be derived from the appropriate particle system. This model is studied through the use of recently developed and benchmarked spectral schemes for the simulation of solutions to stochastic partial differential equations. The results included here demonstrate the effect of adjusting the interparticle interaction on the morphological evolution of the system at the macroscopic level.

REPORT 0506-48: Prediction of mRNA Polyadenylation Sites by Support Vector Machine

¹ Department of Mathematical Sciences, New Jersey Institute of Technology, Newark, NJ 07102

² Department of Biochemistry and Molecular Biology, New Jersey Medical School, University of Medicine and Dentistry of New Jersey, Newark, NJ 07101

Availability: The method is implemented in the program polya_svm, which can be downloaded from http://exon.umdnj.edu/polya_svm.

TECHNICAL REPORTS of the

Center for Applied Mathematics and Statistics