I used a system of ordinary differential equations (ODE) to model a pair of cortical columns, each containing a pool of excitatory and a pool of inhibitory neurons.
I started by considering a single column, with only GABA and NMDA (neurotransmitters) dynamics, so that phase plane methods and simple Hopf bifurcation analysis can be used to find periodic orbits. I then used center-manifold reduction methods on the 3-dimensional AMPA-NMDA-GABA system, confirming the results of the 2-dimensional analysis, and proving the existence of a stable periodic orbit corresponding to a state in which firing rates are modulated periodically.
Besides the computations and proofs done by hand, I used Mathematica and Matlab to solve the systems of ODEs, bring the system into the standard form with a suitable coordinate change (such that the linear part is block diagonal), and compute the center manifold approximation and stability coefficients.