Introduction¶

Most of us are familiar with the ziggagging path of a fly, a behavior has been studied in and out of the lab for many years. Though most of our knowledge comes from animals flying in enclosed environments, in many species we often see a common behavioral motief where straight flight segments are interspersed with sharp turns [mike land]. These sharp turns can be rather stereotyped, and have been termed flight- or body-saccades in analogy with the fast eye movements (oculomotor saccades) performed by vertebrates.

Fortunately, flies will beat their wings when tethered to a small metal pin, a fact that makes it much easier to study their behavior. Even though tethered flight and free flight differ in many significant ways, we can still detect fictive turns in this setup (either by measuring the torqe on the pin or from changes in wing movment). Like free-flight these fictive turns are charactarized by stable epochs broken up by short, large amplitude events that bear the hallmarks of free-flight saccades. Moreover, if we allow the fly to rotate around a single degree of freedom by suspending the tether between two magnets, we will again observe a saccadic turns -- but this time not just from wing motion, but rather from the rotation of the fly within the magnetic field [ref]. A recent paper by Mongeau and Frye combined this 'magno-tether' setup with a panoramic electronic visual display to study how flies use visual motion to trigger saccades, and found evidence for at least two types of behaviors that are distinguished by their dynamics and the aspects of visual motion that drive their initiation. In addition, this study showed how saccades might be integrated with another turning strategy that involves smooth continuous tracking of panoramic motion.

There are several reasons why flies may utilize an episodic vs continuous turning strategies. Many of the proposed mechanisms relate to aspects of sensory processing. For instance, by restricting turns to large, short-duration events the fly might minimize the time in which body rotation dominates the visual motion on their retina. This is a useful feature, since in addition to a system for performing volitional turns, flies must also have many negative-feedback circuits that counteract body rotation to support stability during flight. Restricting the duration of volitional turns would help a turn-promoting saccade system and a turn-suppressing stabilization system coexist within a single fly brain, since they need not be active at the same time. Consistent with this idea, there is physiological evidence for a circuit that suppresses visual feedback during saccades. Moreover, there is even evidence for an efferent system that could modulate the sensitivity of mechanosensory feedback during saccades, but how this embodied control circut interfaces with saccade production is still being worked out.

In addition to these sensory explanations for the existence of saccades, the behavior may be adaptive from the perspective of the motor system as well. The wing is controlled by a very sparse set of twelve muscles, each innervated by only one motor neuron. In my most recent paper, I used a fixed-tether setup combined with calcium imaging to study how flies recruit these flight muscles to control wing motion, and found that they could be functionally separated into two groups: one group of smaller muscles that are persistently active, and appear to be modulated during slow, smooth changes in wing motion, and a second group of large muscles that are only sporadically active during saccades. This stratification suggests a mechanism where flight saccades might be generated as an emergent property of the muscle recruitment thresholds. Additionally, because of the dynamics of calcium cycling in muscles, this mechanism may provide an explanation for the characteristic time course of tethered flight saccades with their rapid onset and slower exponential decay. In the figure below, you can see how wing kinematics during saccades tracks the activity of the underlying muscles, i1 and iii3. It is important to take this result with caution however, because sensor we use to monitor muscle activity is a little slow, and may not accurately follow the underlying activity with perfect fidelity. Nevertheless, if we simply fit the wing amplitude as the sum of the scaled activity signals we can explain anywhere from 70-80% of the variation in wing kinematics at this temporal resolution.

To clarify the quantitative implications of this mechansism, I generated a simple computational (nlp) model of the flight motor system built around a bank of motor neurons driven by noisy inputs. The characteristic feature of the model comes down to differences in the input-output functions of the motor neurons: the sporadically active muscles (which we refer to as phasic) have a steeper rightward shifted threshold, whereas the persistently active muscles (which we refer to as tonic) have a shallower leftward shifted threshold. These threshold differences mean that the phasic muscles are only recruited by inputs that come from the tails of the noise distribution and thus only fire brief bursts of spikes. In contrast, the tonic muscles thresholds are centered with respect to the noisy input and therefore respond by making rapid transitions between the floor and cieling of the input output function.

There are a number of free parameters in this model, which I have adjusted by hand with the intention of showing that this simple mechanism is:

Sufficient to recapitulate many of the phenomena that we observe in our experiments.
Built from well established neuronal computations.
Incorporates what is known about the sensory and descending inputs into the motor system.

It would be exciting in future work to try and fit the parameters of the model to our calcium imaging data, or perhaps directly measure membrane properties and synaptic inputs from electrical recordings of the motor neurons.

In [31]:

TONIC_BCK = 1.0#250 # background input from the tonic channels
PHASIC_BCK = 0# 20 # background input from the phasic channels

INT_GAIN = 60 #gain on the inputs from the visual motion pathway
SAC_GAIN = 12 #gain on the noisy input driving saccades
#INT_GAIN = 0 #set to zero to study the noise process in isoloation 

SK_TON = 50 #on time-constant of the visual integrator kernel
SK_TOFF = 5000 #off time-constant of the visual integrator kernel

MK_TON = 10 #on time-constant of muscle calcium dynamics
MK_TOFF = 200 #off time-constant of muscle calcium dynamics

In [32]:

# setup the environment
%matplotlib inline
%config InlineBackend.figure_format = 'png' 
import pylab as plb
import numpy as np
import scipy as sp
import figurefirst as fifi
#import plotfuncs as plf
from numpy import random
plb.rcParams['font.size'] =  9
plb.rcParams['pdf.fonttype'] = 42
plb.rcParams['image.interpolation'] = 'nearest'
plb.rcParams['image.aspect'] = 'auto'
plb.rcParams['image.cmap'] = 'BrBG'
demo_mode = True
if demo_mode:
    random.seed(2)

In [33]:

# Timebase
xi = np.linspace(-20,20,40000) #times
sf = 40000/40.0 #sampling frequency
step = ((xi>0) & (xi<3)).astype(float) #motion stimulus
rng = ((xi>-3) & (xi<6)) # the plot range

In [34]:

def make_kernel(tauon,tauoff):
    kx = np.linspace(0.1,20000,10000)
    kon = lambda x:np.exp(((-1*tauon)/(x)))
    koff = lambda x:np.exp((-1*x)/tauoff)
    k = (kon(kx)*koff(kx)/np.max(kon(kx)*koff(kx)))
    k = np.hstack((np.zeros_like(k),k)) #acausal
    return k

sk = make_kernel(SK_TON,SK_TOFF) #slow kernal (optimotor integrator see below)
mk = make_kernel(MK_TON,MK_TOFF) #muscle kernel -- calcium and GcAMP

plb.plot(np.linspace(-20,20,20000),mk,'g')
plb.plot(np.linspace(-20,20,20000),sk,'r')
plb.gca().set_xlabel('time (s)')
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])

Static nonlinearity¶

The phasic muscles have a sharper and higher activation threshold than the tonic muscles. I model this as differences in a static non-linearity that defines the input-output threshold of each motor neuron. Note that this could arise via differences in neuronal excitability, the strengths of the synaptic inputs, or both. It is known that some if not all of these motor neurons receive a strong stroke-by-stroke synaptic impulse from mechanosensors on the wing and haltere (a modified hindwing with a sensory function). For this reason, one could imagine that the form of the input-output functions are shaped by the strength and precision of these mechanosensory inputs. Additionally, since the fly wingstroke cycle is driven by an indirect myogenic oscillator at 200 Hz these motor neurons can really only fire one spike per wingstroke, therefore I can set the maximum spike frequency to wingstroke frequency.

In [35]:

#muscle firing rate can range from 0 to 200 Hz
sigmoid_t = lambda x: 200/(1+np.exp(-0.03*(x - 150))) #tonic squashing function
sigmoid_p = lambda x: 200/(1+np.exp(-0.08*(x - 400))) #phasic squashing function
inp = np.linspace(0,600,100) #plot example sigmoid
plb.plot(inp,sigmoid_t(inp))
plb.plot(inp,sigmoid_p(inp))
plb.gca().set_ylabel('output rate')
plb.gca().set_xlabel('input rate')
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])

Noise source¶

I am modeling a common source of noise to the entire system. For simplicity, I am using the sum of a positive and negative poisson process, so that the net input has a mean of zero. The statistics of this noise process translate into statistics of saccade behaviors: the heavier the tail the more infrequent are saccades. The underlying process will also determine the inter-saccade interval - a Poisson process predicts an exponential distribution, whereas a power law will yield a levy flight.

In [36]:

# create a common saccade trigger signal as the sum
# of a positive and negative poisson process
BK_RATE = 60000.0 # background rate
#MOTION_RATE = 5.0 # rate during the motion epoch
n_common = np.random.poisson((np.zeros_like(step)+BK_RATE)/sf)
n_common -= np.random.poisson((np.zeros_like(step)+BK_RATE)/sf)
plb.axvspan(0,3,color = 'k',alpha = 0.3)
plb.plot(xi[rng],n_common[rng])
plb.gca().set_ylabel('rightward rate             leftward rate')
plb.gca().set_xlabel('time (s)')
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])

Optimotor signal¶

We know that flies combine a smooth response to wide-field motion with the production of saccades, therefore, the motor neurons must also be able to incorporate visual inputs. Experiments have measured the time-course of the behavioral response to open-loop motion and found that the response has the characteristics of a leaky integrator. To model this, I take the open loop motion stimulus and run it through a low-pass filter, the kernel is shown in red above, and is very slow in comparison to the 'muscle kernel'. The signal from the optimotor pathway is a signed in a way that depends on the direction of motion, in this case it is positive for leftward motion. The synaptic weights are such that a positive signal will excite the motor neurons that increase right wing amplitude and inhibit the motor neurons that decrease left wing amplitude.

In [37]:

# create common integrated optomotor signal by convovling the motion step with the 
# slow kernel
integrate = np.convolve(sk,step,mode = 'same')
integrate /= max(integrate)
plb.axvspan(0,3,color = 'k',alpha = 0.3)
plb.plot(xi[rng],integrate[rng])
#plb.plot(xi[rng],-1*integrate[rng])
plb.gca().set_ylabel('leftward drive')
plb.gca().set_xlabel('time (s)')
plb.gca().set_xticks(np.arange(-3,6))
plb.gca().set_xticklabels(np.arange(-3,6))
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])

Motoneuron processing¶

The motor unit activation signal is a modeled by summing the descending signals, passing them through a static non-linearity that determines the time-varying motor neuron firing rate. This firing rate is then filtered by the muscle calcium dynamics.

The muscle array is represented by four types of motor units for each wing: a phasic and tonic unit for the muscles associated with increases in stroke amplitude (blue), and a phasic and tonic unit for the muscles associated with a decrease in stroke amplitude (green). The ampup and ampdown muscles differ in that the polarity of the inputs are reversed. The input rates for the muscles attached to the right wing are shown below.

In [38]:

rwing_phasic_ampup_f = sigmoid_p(integrate*INT_GAIN + (n_common+ PHASIC_BCK)*SAC_GAIN )
rwing_phasic_ampdown_f = sigmoid_p(-1*integrate*INT_GAIN - (n_common+ PHASIC_BCK)*SAC_GAIN )
lwing_phasic_ampup_f = sigmoid_p(-1*integrate*INT_GAIN - (n_common+ PHASIC_BCK)*SAC_GAIN )
lwing_phasic_ampdown_f = sigmoid_p(integrate*INT_GAIN + (n_common+ PHASIC_BCK)*SAC_GAIN )

rwing_tonic_ampup_f = sigmoid_t(integrate*INT_GAIN + (n_common+ TONIC_BCK)*SAC_GAIN  )
rwing_tonic_ampdown_f = sigmoid_t(-1*integrate*INT_GAIN - (n_common+ TONIC_BCK)*SAC_GAIN  )
lwing_tonic_ampup_f = sigmoid_t(-1*integrate*INT_GAIN - (n_common+ TONIC_BCK)*SAC_GAIN  )
lwing_tonic_ampdown_f = sigmoid_t(integrate*INT_GAIN + (n_common + TONIC_BCK)*SAC_GAIN )

plb.subplot(2,1,1)
plb.plot(xi[rng],rwing_phasic_ampup_f[rng],color = 'b')
plb.plot(xi[rng],rwing_phasic_ampdown_f[rng],color = 'g')
plb.axvspan(0,3,color = 'k',alpha = 0.2)
plb.gca().set_ylabel('phasic poisson rate')
fifi.mpl_functions.adjust_spines(plb.gca(),['left'])

plb.subplot(2,1,2)
plb.plot(xi[rng],rwing_tonic_ampup_f[rng],color = 'b')
plb.plot(xi[rng],rwing_tonic_ampdown_f[rng],color = 'g')
plb.axvspan(0,3,color = 'k',alpha = 0.2)
plb.gca().set_ybound(0,250)
plb.gca().set_ylabel('tonic poisson rate')
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])
lbl = plb.gca().set_xlabel('time (s)')

In [39]:

# to create the calcium signals we pass the muscle firing rates through the 'muscle kernel'
ca_sigs = dict()
ca_sigs[('right_wing','phasic','amp_up')] = np.convolve(mk,rwing_phasic_ampup_f,mode = 'same')
ca_sigs[('right_wing','phasic','amp_down')] = np.convolve(mk,rwing_phasic_ampdown_f,mode = 'same')
ca_sigs[('left_wing','phasic','amp_up')] = np.convolve(mk,lwing_phasic_ampup_f,mode = 'same')
ca_sigs[('left_wing','phasic','amp_down')] = np.convolve(mk,lwing_phasic_ampdown_f,mode = 'same')

ca_sigs[('right_wing','tonic','amp_up')] = np.convolve(mk,rwing_tonic_ampup_f,mode = 'same')
ca_sigs[('right_wing','tonic','amp_down')] = np.convolve(mk,rwing_tonic_ampdown_f,mode = 'same')
ca_sigs[('left_wing','tonic','amp_up')] = np.convolve(mk,lwing_tonic_ampup_f,mode = 'same')
ca_sigs[('left_wing','tonic','amp_down')] = np.convolve(mk,lwing_tonic_ampdown_f,mode = 'same')

for sig in ca_sigs.values():
    sig -= np.mean(sig[(xi<0)&(xi>-3)])
    sig /= np.std(sig[(xi<0)&(xi>-3)])
    
plb.subplot(2,1,1)

plb.plot(xi[rng],ca_sigs[('right_wing','phasic','amp_up')][rng])
plb.plot(xi[rng],ca_sigs[('right_wing','phasic','amp_down')][rng])
plb.axvspan(0,3,color = 'k',alpha = 0.2)
fifi.mpl_functions.adjust_spines(plb.gca(),['left'])
plb.gca().set_ylabel('right wing activation')

plb.subplot(2,1,2)
plb.plot(xi[rng],ca_sigs[('right_wing','tonic','amp_up')][rng])
plb.plot(xi[rng],ca_sigs[('right_wing','tonic','amp_down')][rng])
plb.axvspan(0,3,color = 'k',alpha = 0.2)
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])
plb.gca().set_ylabel('left wing activation')
plb.gca().set_xlabel('time (s)')
lbl = plb.gcf().suptitle('comparison of phasic muscle activity on left vs right wing')

Wing kinematics¶

To model the wing kinematics and associated body torqe, I add up the tonic and phasic on signals on the right and left wing to get the right and left wingstroke amplitude respectively. The turning torque is modeled by subtracting the right amplitude from the left amplitude. Note that the left and right wingstroke amplitude are highly anticorelated with one another, a result of the fact that the same noise process is driving the left and right wing.

In [40]:

rw = 0.5*ca_sigs[('right_wing','tonic','amp_up')] + 1*ca_sigs[('right_wing','phasic','amp_up')] - \
            0.5*ca_sigs[('right_wing','tonic','amp_down')]-1*ca_sigs[('right_wing','phasic','amp_down')] 
    
lw = 0.5*ca_sigs[('left_wing','tonic','amp_up')] + 1*ca_sigs[('left_wing','phasic','amp_up')] - \
        0.5*ca_sigs[('left_wing','tonic','amp_down')]-1*ca_sigs[('left_wing','phasic','amp_down')]
lmr = lw-rw

In [41]:

plb.plot(xi[rng],lw[rng])
plb.plot(xi[rng],rw[rng])
plb.plot(xi[rng],lw[rng]-rw[rng])
plb.axvspan(0,3,color = 'k',alpha = 0.2)
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])
plb.gca().set_yticklabels([])
plb.gca().set_xlabel('time (s)')
lbl = plb.gca().set_ylabel('wing amplitude (AU)')

Summary figures¶

To summarize the output of the model for the given parameters I will use the figurefirst layout document 'model_detal_layout.svg'. I will also generate another figure using the 'quant_outline_layout.svg' that compares the details of the model a tonic vs phasic model.

In [12]:

!wget 'https://storage.googleapis.com/lindsay-labbook-public-images/model_detail_layout.svg'

--2018-03-22 22:14:44--  https://storage.googleapis.com/lindsay-labbook-public-images/model_detail_layout.svg
Resolving storage.googleapis.com (storage.googleapis.com)... 74.125.124.128, 2607:f8b0:4001:c12::80
Connecting to storage.googleapis.com (storage.googleapis.com)|74.125.124.128|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 47358 (46K) [image/svg+xml]
Saving to: ‘model_detail_layout.svg.3’

model_detail_layout 100%[===================>]  46.25K  --.-KB/s    in 0.001s  

2018-03-22 22:14:44 (84.2 MB/s) - ‘model_detail_layout.svg.3’ saved [47358/47358]

In [13]:

!wget 'https://storage.googleapis.com/lindsay-labbook-public-images/quant_outline_layout.svg'

--2018-03-22 22:14:44--  https://storage.googleapis.com/lindsay-labbook-public-images/quant_outline_layout.svg
Resolving storage.googleapis.com (storage.googleapis.com)... 74.125.124.128, 2607:f8b0:4001:c12::80
Connecting to storage.googleapis.com (storage.googleapis.com)|74.125.124.128|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 443485 (433K) [image/svg+xml]
Saving to: ‘quant_outline_layout.svg.3’

quant_outline_layou 100%[===================>] 433.09K  --.-KB/s    in 0.003s  

2018-03-22 22:14:44 (129 MB/s) - ‘quant_outline_layout.svg.3’ saved [443485/443485]

In [42]:

reload(fifi)
from matplotlib import colors
layout = fifi.FigureLayout('quant_outline_layout.svg')
layout.make_mplfigures()

rng = ((xi>-3) & (xi<6))

tmax = np.max(np.hstack([ca_sigs[('left_wing','tonic','amp_up')][rng], 
                         ca_sigs[('left_wing','tonic','amp_down')][rng],
                         ca_sigs[('right_wing','tonic','amp_up')][rng],
                         ca_sigs[('right_wing','tonic','amp_down')][rng]]))

pmax = np.max(np.hstack([ca_sigs[('left_wing','phasic','amp_up')][rng], 
                         ca_sigs[('left_wing','phasic','amp_down')][rng],
                         ca_sigs[('right_wing','phasic','amp_up')][rng],
                         ca_sigs[('right_wing','phasic','amp_down')][rng]]))


for wing_key in ['left_wing','right_wing']:
    for dyn_key in ['tonic','phasic']:
        for amp_key in ['amp_up','amp_down']:
            ax = layout.axes[('model',wing_key,dyn_key,amp_key)]
            ax.plot(xi[rng],ca_sigs[(wing_key,dyn_key,amp_key)][rng],clip_on = False,color = 'k')
            if dyn_key == 'tonic':
                ax.set_ybound(-1*tmax,tmax)
            else:
                ax.set_ybound(0,pmax)

### wing_kine
mx = np.max(lmr[rng]);mn = np.min(lmr[rng])
rn = mx-mn;
mx = rn/2.0
mn = -rn/2.0
ax = layout.axes_groups['none']['model']['YL']['left_wing']
ax.plot(xi[rng],lw[rng],clip_on = False,color = 'k')
ax.set_ybound(mn,mx)

ax = layout.axes_groups['none']['model']['YL']['right_wing']
ax.plot(xi[rng],rw[rng],clip_on = False,color = 'k')
ax.set_ybound(mn,mx)

ax = layout.axes_groups['none']['model']['YL']['lmr']
ax.plot(xi[rng],(lw-rw)[rng],clip_on = False,color = 'k')
ax.set_ybound(mn,mx)

for key,ax in layout.axes.items():
    if (key[0] == 'model') & ~(key[1] =='squashing'):
        plb.sca(ax['axis'])
        ax.set_xbound(-3,6)
        #ax.set_ybound(-1,1)
        if INT_GAIN:
            plb.axvspan(0,3,alpha = 0.1,color = 'k')
    fifi.mpl_functions.kill_spines(ax)
        #plf.kill_spines()
        
        
### sigmoids
for axkey,ax in layout.axes_groups['none']['model']['squashing']['left_squash'].items():
    if 't' in axkey:
        if 'up' in axkey:
            shift = -45
        if 'down' in axkey:
            shift = 45
        sigmoid = lambda x: 300/(1+np.exp(-0.02*(x - 150))) #squashing function
        inp_ = np.linspace(0,300,100) #plot example sigmoid
        ax.plot(inp_,sigmoid(inp_),color = 'k',clip_on = False)
        ax.axhline(sigmoid(150+shift),ls = (0,(0.5,0.5)),color = 'k')
        ax.axvline(150+shift,ls = (0,(0.5,0.5)),color = 'k')
    else:# the plot range
        if 'up' in axkey:
            shift = -15
        if 'down' in axkey:
            shift = -5
        sigmoid = lambda x: 300/(1+np.exp(-0.1*(x - 150))) #squashing function
        inp_ = np.linspace(0+shift,300+shift,100)
        ax.plot(inp_,sigmoid(inp_),color = 'k',clip_on = False)
        ax.axhline(sigmoid(150+shift),ls = (0,(0.5,0.5)),color = 'k')
        ax.axvline(150+shift,ls = (0,(0.5,0.5)),color = 'k')

        
for axkey,ax in layout.axes_groups['none']['model']['squashing']['right_squash'].items():
    if 't' in axkey:
        if 'up' in axkey:
            shift = 45
        if 'down' in axkey:
            shift = -45
        sigmoid = lambda x: 300/(1+np.exp(-0.02*(x - 150))) #squashing function
        inp_ = np.linspace(0,300,100) #plot example sigmoid
        ax.plot(inp_,sigmoid(inp_),color = 'k',clip_on = False)
        ax.axhline(sigmoid(150+shift),ls = (0,(0.5,0.5)),color = 'k')
        ax.axvline(150+shift,ls = (0,(0.5,0.5)),color = 'k')
    else:
        if 'up' in axkey:
            shift = -5
        if 'down' in axkey:
            shift = -15
        sigmoid = lambda x: 300/(1+np.exp(-0.1*(x - 150))) #squashing function
        inp_ = np.linspace(0+shift,300+shift,100)
        ax.plot(inp_,sigmoid(inp_),color = 'k',clip_on = False)
        ax.axhline(sigmoid(150+shift),ls = (0,(0.5,0.5)),color = 'k')
        ax.axvline(150+shift,ls = (0,(0.5,0.5)),color = 'k')
if INT_GAIN:
    layout.axes_groups['none']['model']['int_command'].plot(xi[rng],integrate[rng],color = 'k')
layout.axes_groups['none']['model']['sac_command'].plot(xi[rng],n_common[rng],color = 'k')


plb.close('all')
layout.insert_figures()
#layout.apply_svg_attrs()
layout.set_layer_visibility('Layer 1',False)
layout.write_svg('quant_outline.svg')
from IPython.display import display,SVG
display(SVG('quant_outline.svg'))

In [28]:

mpl.rcParams.keys()#['axes.autolimit_mode']

Out[28]:

[u'agg.path.chunksize',
 u'animation.avconv_args',
 u'animation.avconv_path',
 u'animation.bitrate',
 u'animation.codec',
 u'animation.convert_args',
 u'animation.convert_path',
 u'animation.ffmpeg_args',
 u'animation.ffmpeg_path',
 u'animation.frame_format',
 u'animation.html',
 u'animation.mencoder_args',
 u'animation.mencoder_path',
 u'animation.writer',
 u'axes.axisbelow',
 u'axes.edgecolor',
 u'axes.facecolor',
 u'axes.formatter.limits',
 u'axes.formatter.use_locale',
 u'axes.formatter.use_mathtext',
 u'axes.formatter.useoffset',
 u'axes.grid',
 u'axes.grid.axis',
 u'axes.grid.which',
 u'axes.hold',
 u'axes.labelcolor',
 u'axes.labelpad',
 u'axes.labelsize',
 u'axes.labelweight',
 u'axes.linewidth',
 u'axes.prop_cycle',
 u'axes.spines.bottom',
 u'axes.spines.left',
 u'axes.spines.right',
 u'axes.spines.top',
 u'axes.titlesize',
 u'axes.titleweight',
 u'axes.unicode_minus',
 u'axes.xmargin',
 u'axes.ymargin',
 u'axes3d.grid',
 u'backend',
 u'backend.qt4',
 u'backend.qt5',
 u'backend_fallback',
 u'boxplot.bootstrap',
 u'boxplot.boxprops.color',
 u'boxplot.boxprops.linestyle',
 u'boxplot.boxprops.linewidth',
 u'boxplot.capprops.color',
 u'boxplot.capprops.linestyle',
 u'boxplot.capprops.linewidth',
 u'boxplot.flierprops.color',
 u'boxplot.flierprops.linestyle',
 u'boxplot.flierprops.linewidth',
 u'boxplot.flierprops.marker',
 u'boxplot.flierprops.markeredgecolor',
 u'boxplot.flierprops.markerfacecolor',
 u'boxplot.flierprops.markersize',
 u'boxplot.meanline',
 u'boxplot.meanprops.color',
 u'boxplot.meanprops.linestyle',
 u'boxplot.meanprops.linewidth',
 u'boxplot.medianprops.color',
 u'boxplot.medianprops.linestyle',
 u'boxplot.medianprops.linewidth',
 u'boxplot.notch',
 u'boxplot.patchartist',
 u'boxplot.showbox',
 u'boxplot.showcaps',
 u'boxplot.showfliers',
 u'boxplot.showmeans',
 u'boxplot.vertical',
 u'boxplot.whiskerprops.color',
 u'boxplot.whiskerprops.linestyle',
 u'boxplot.whiskerprops.linewidth',
 u'boxplot.whiskers',
 u'contour.corner_mask',
 u'contour.negative_linestyle',
 u'datapath',
 u'docstring.hardcopy',
 u'errorbar.capsize',
 u'examples.directory',
 u'figure.autolayout',
 u'figure.dpi',
 u'figure.edgecolor',
 u'figure.facecolor',
 u'figure.figsize',
 u'figure.frameon',
 u'figure.max_open_warning',
 u'figure.subplot.bottom',
 u'figure.subplot.hspace',
 u'figure.subplot.left',
 u'figure.subplot.right',
 u'figure.subplot.top',
 u'figure.subplot.wspace',
 u'figure.titlesize',
 u'figure.titleweight',
 u'font.cursive',
 u'font.family',
 u'font.fantasy',
 u'font.monospace',
 u'font.sans-serif',
 u'font.serif',
 u'font.size',
 u'font.stretch',
 u'font.style',
 u'font.variant',
 u'font.weight',
 u'grid.alpha',
 u'grid.color',
 u'grid.linestyle',
 u'grid.linewidth',
 u'image.aspect',
 u'image.cmap',
 u'image.composite_image',
 u'image.interpolation',
 u'image.lut',
 u'image.origin',
 u'image.resample',
 u'interactive',
 u'keymap.all_axes',
 u'keymap.back',
 u'keymap.forward',
 u'keymap.fullscreen',
 u'keymap.grid',
 u'keymap.home',
 u'keymap.pan',
 u'keymap.quit',
 u'keymap.save',
 u'keymap.xscale',
 u'keymap.yscale',
 u'keymap.zoom',
 u'legend.borderaxespad',
 u'legend.borderpad',
 u'legend.columnspacing',
 u'legend.edgecolor',
 u'legend.facecolor',
 u'legend.fancybox',
 u'legend.fontsize',
 u'legend.framealpha',
 u'legend.frameon',
 u'legend.handleheight',
 u'legend.handlelength',
 u'legend.handletextpad',
 u'legend.isaxes',
 u'legend.labelspacing',
 u'legend.loc',
 u'legend.markerscale',
 u'legend.numpoints',
 u'legend.scatterpoints',
 u'legend.shadow',
 u'lines.antialiased',
 u'lines.color',
 u'lines.dash_capstyle',
 u'lines.dash_joinstyle',
 u'lines.linestyle',
 u'lines.linewidth',
 u'lines.marker',
 u'lines.markeredgewidth',
 u'lines.markersize',
 u'lines.solid_capstyle',
 u'lines.solid_joinstyle',
 u'markers.fillstyle',
 u'mathtext.bf',
 u'mathtext.cal',
 u'mathtext.default',
 u'mathtext.fallback_to_cm',
 u'mathtext.fontset',
 u'mathtext.it',
 u'mathtext.rm',
 u'mathtext.sf',
 u'mathtext.tt',
 u'nbagg.transparent',
 u'patch.antialiased',
 u'patch.edgecolor',
 u'patch.facecolor',
 u'patch.linewidth',
 u'path.effects',
 u'path.simplify',
 u'path.simplify_threshold',
 u'path.sketch',
 u'path.snap',
 u'pdf.compression',
 u'pdf.fonttype',
 u'pdf.inheritcolor',
 u'pdf.use14corefonts',
 u'pgf.debug',
 u'pgf.preamble',
 u'pgf.rcfonts',
 u'pgf.texsystem',
 u'plugins.directory',
 u'polaraxes.grid',
 u'ps.distiller.res',
 u'ps.fonttype',
 u'ps.papersize',
 u'ps.useafm',
 u'ps.usedistiller',
 u'savefig.bbox',
 u'savefig.directory',
 u'savefig.dpi',
 u'savefig.edgecolor',
 u'savefig.facecolor',
 u'savefig.format',
 u'savefig.frameon',
 u'savefig.jpeg_quality',
 u'savefig.orientation',
 u'savefig.pad_inches',
 u'savefig.transparent',
 u'svg.fonttype',
 u'svg.image_inline',
 u'svg.image_noscale',
 u'text.antialiased',
 u'text.color',
 u'text.dvipnghack',
 u'text.hinting',
 u'text.hinting_factor',
 u'text.latex.preamble',
 u'text.latex.preview',
 u'text.latex.unicode',
 u'text.usetex',
 u'timezone',
 u'tk.window_focus',
 u'toolbar',
 u'verbose.fileo',
 u'verbose.level',
 u'webagg.open_in_browser',
 u'webagg.port',
 u'webagg.port_retries',
 u'xtick.color',
 u'xtick.direction',
 u'xtick.labelsize',
 u'xtick.major.pad',
 u'xtick.major.size',
 u'xtick.major.width',
 u'xtick.minor.pad',
 u'xtick.minor.size',
 u'xtick.minor.visible',
 u'xtick.minor.width',
 u'ytick.color',
 u'ytick.direction',
 u'ytick.labelsize',
 u'ytick.major.pad',
 u'ytick.major.size',
 u'ytick.major.width',
 u'ytick.minor.pad',
 u'ytick.minor.size',
 u'ytick.minor.visible',
 u'ytick.minor.width']

In [29]:

In [15]:

rng = ((xi>-3) & (xi<6))
i_rng = ((xi>0) & (xi<0.5))

ca_sigs = dict()
ca_sigs[('right_wing','phasic','amp_up')] = np.convolve(mk,rwing_phasic_ampup_f,mode = 'same')
ca_sigs[('right_wing','phasic','amp_down')] = np.convolve(mk,rwing_phasic_ampdown_f,mode = 'same')
ca_sigs[('left_wing','phasic','amp_up')] = np.convolve(mk,lwing_phasic_ampup_f,mode = 'same')
ca_sigs[('left_wing','phasic','amp_down')] = np.convolve(mk,lwing_phasic_ampdown_f,mode = 'same')

ca_sigs[('right_wing','tonic','amp_up')] = np.convolve(mk,rwing_tonic_ampup_f,mode = 'same')
ca_sigs[('right_wing','tonic','amp_down')] = np.convolve(mk,rwing_tonic_ampdown_f,mode = 'same')
ca_sigs[('left_wing','tonic','amp_up')] = np.convolve(mk,lwing_tonic_ampup_f,mode = 'same')
ca_sigs[('left_wing','tonic','amp_down')] = np.convolve(mk,lwing_tonic_ampdown_f,mode = 'same')

#rng = ((xi>-1) & (xi<3))
layout = fifi.FigureLayout('model_detail_layout.svg',make_mplfigures = True)
p_color = layout.pathspecs['phasic'].mplkwargs()['edgecolor']
t_color = layout.pathspecs['tonic'].mplkwargs()['edgecolor']

layout.axes['poisson'].plot(xi[rng],n_common[rng],color = 'k')

layout.axes['thresholds_t'].plot(inp/SAC_GAIN,sigmoid_t(inp),
                                 clip_on = False,color = t_color)
layout.axes['thresholds_p'].plot(inp/SAC_GAIN,sigmoid_p(inp),
                                 clip_on = False,color = p_color)

layout.axes['tonic_inset'].plot(xi[i_rng],
                                rwing_tonic_ampup_f[i_rng],
                               color = t_color)

layout.axes['tonic_rate'].plot(xi[rng],
                               rwing_tonic_ampup_f[rng],
                               color = t_color)
layout.axes['phasic_rate'].plot(xi[rng],
                                rwing_phasic_ampup_f[rng],
                                color = p_color)


layout.axes['calcium_kernel_p'].plot(np.linspace(-20,20,20000),mk,color = p_color)

layout.axes['calcium_kernel_p'].set_xbound(0,1.0)
layout.axes['calcium_kernel_t'].plot(np.linspace(-20,20,20000),mk,color = t_color)
layout.axes['calcium_kernel_t'].set_xbound(0,1.0)
layout.axes['calcium_kernel_t'].set_yticks([0,0.5,1.0])
layout.axes['calcium_kernel_p'].set_yticks([0,0.5,1.0])
layout.axes['calcium_kernel_p'].set_xticks([0,0.5,1.0])

layout.axes['tonic_calcium'].plot(xi[rng],
                                  ca_sigs[('left_wing','tonic','amp_up')][rng],
                                  color = t_color)
layout.axes['phasic_calcium'].plot(xi[rng],
                                   ca_sigs[('left_wing','phasic','amp_up')][rng],
                                   color = p_color)

fifi.mpl_functions.set_spines(layout)
layout.save('model_detail.svg')
plb.close('all')
display(SVG('model_detail.svg'))
#plb.plot(inp,sigmoid_p(inp))

Summary¶

This model is by no means complete, but I do think it highlights a few reasons why it is important to consider biomechanics when trying to understand how the nervous system controls behavior. First, it is worth recognizing that the sparsity of the steering array in Drosophila may play into the phenomena we recognize as a flight saccade; with only a limited set of control elements, the fly may be forced into recruiting muscles that have a relatively large effect on wing motion (the phasic muscles) from time-to-time, and thus generate a stereotyped impulsive turn. Second, once the fly recruits these muscles -- even if for a short period -- calcium may remain elevated for a number of wingstrokes after the cell has stopped spiking. In this way, the fly would be 'stuck' with partially activated phasic muscles until calcium levels can decay to baseline. For a tethered fly, this could provide an explanation for the duration and time-course of the observed events. As can be seen in left-right trace in the first figure, the tethered-flight saccades decay to baseline over the course of a couple hundred milliseconds. In contrast, a free-flight saccade is completed much faster, 2-3 wingstrokes or around 10-20ms. The reasons for this discrepancy in duration are still not fully understood; however, magnotether experiments indicate that a fast mechanosensory feedback pathway that detects body rotation may serve to accelerate the termination of saccades when the fly is not fixed firmly in place.

The model I propose here predicts that the motor program the fly uses during free flight saccades would need to take into account muscle calcium dynamics to successfully end the turn. To be explicit, the fly could not simply wait for the muscles it used to initiate the saccade to turn off, but would instead rely on faster mechanisms such as the recruitment of antagonistic muscles or subtle adjustments in the timing of spikes within the stroke cycle [ref]. A similar situation may exist within the asynchronous motor system that powers the wingstroke, where calcium dynamics are likely even slower than in the steering muscles.

To test these implications, a key direction for future work would be to identify the correlations and sequences within the motor array at higher temporal resolution. Clearly, the relative importance of muscle dynamics is highly dependent on how fast these steering muscles can regulate their activity. In the scenario I model, muscle calcium cycling is relatively slow, such that a phasic muscle would driven by a short burst of spikes at the onset of a saccade and the decay in the induced wing kinematics end up tracking calcium (a process we think the nervous system has little-to-no control over). In contrast, if muscle calcium is relatively fast, we would expect the hidden, underlying spike rate to decay as the saccade winds down . As a final possibility, one that is compatible with both slow and fast muscle cycling, the saccade time course may be build from a staggered temporal sequence of muscle recruitment across the population. This staggering, might occur within the muscles of a given wing, or accross the left and right wings.

So how fast is calcium dynamics in muscle? Previous in-vivo work loop experiments in blowflies has shown that when stimulated with a 150Hz pulse train, the b1 muscle (a tonic muscle) doesn't completely reach tetanus, an important findings since the remaining ripple of force provides a fast mechanism for stroke-by-stroke regulation of muscle work. It is unclear however, if this is would hold true for fruit-flies which flap their wings at a higher frequency - closer to 200Hz. Compared to vertebrate muscle, insect muscle is ultrastructurally very diverse, raising the possibility that calcium dynamics may vary on a muscle-by-muscle basis. Indeed, the dynamics of a force transient in cicada tymbal muscles are highly correlated with cellular volume occupied by the t-tubule and sarcoplasmic reticulum. This type of ultrastructural analysis has yet to be done for fruit fly flight muscles; however, because flight requires specialization over a wide range of performance variables, a high level of diversity is expected. Looking to the future, the expanding toolkit of Drosophila genetics may give us hope in working out how the how the nervous system and the muscles are co-adapted to produce the deft in-flight turns of this little animal.

In [14]:

layout = fifi.FigureLayout('competing_dynamics_layout.svg',make_mplfigures=True)

In [16]:

from IPython.display import display,SVG
display(SVG('competing_dynamics_layout.svg'))l

In [24]:

SAC_TON = 50 #on time-constant of descending command
SAC_TOFF = 500 #off time-constant of descending command

MKF_TON = 10 #on time-constant of muscle calcium dynamics
MKF_TOFF = 50 #off time-constant of muscle calcium dynamics

MKS_TON = 10 #on time-constant of muscle calcium dynamics
MKS_TOFF = 500 #off time-constant of muscle calcium dynamics

sac = make_kernel(SAC_TON,SAC_TOFF) #slow kernal (optimotor integrator see below)
mk_fast = make_kernel(MKF_TON,MKF_TOFF) #muscle kernel -- calcium and GcAMP
mk_slow = make_kernel(MKS_TON,MKS_TOFF) #muscle kernel -- calcium and GcAMP


plb.plot(np.linspace(-20,20,20000),sac,'g')
plb.plot(np.linspace(-20,20,20000),mk_fast,'r')
plb.plot(np.linspace(-20,20,20000),mk_slow,'b')
plb.gca().set_xlabel('time (s)')
fifi.mpl_functions.adjust_spines(plb.gca(),['left','bottom'])

In [ ]:

Generative model of fast flight turns using motor recrutment thresholds.