User Tools

Site Tools


:!: There will be NO class meeting on the 20th, we will discuss this module and the next on the 27th!

Spike train analysis: firing rate, interspike interval distributions, auto- and cross-correlation


  • Understand the need for statistical measures to describe patterns of spikes
  • Gain familiarity with basic descriptors of spike trains: firing rate, interspike interval distribution, autocorrelation function, cross-correlation function
  • Generate artificial Poisson spike trains and compare these to real data
  • Develop intuitions about how Poisson spike trains can be used to test interesting hypotheses about neural coding



  • Explore neural coding in the rat hippocampus by comparing cross-correlations to synthetic data


Action potentials, or spikes, are the main unit of information processing in the brain. Spikes enable fast yet reliable electrical communication over more than the microscopic distances to which passive conduction is limited. Although there are exceptions, spikes are generally the only aspect of neural activity that is transmitted rapidly to postsynaptic targets. At the most general level, sensory input to the brain can be thought of as changing the patterns of spikes in associated sensory nuclei and brain regions; in turn, motor output is ultimately governed by patterns of spikes in motor areas. Thus, any adaptive behavior must ultimately result from the appropriate manipulation of spike patterns. As a result, the recording and interpretation of spike data is a central tool in neuroscience.

This module introduces the most commonly used basic tools for describing spike trains. We focus first on describing spike trains from a single neuron, with a first step towards pairwise analysis at the end.

Estimating firing rate

For most purposes, a spike can be considered as a unitary, punctate event that occurs at a specific time. Of course, in reality, a spike actually has a certain brief lifetime, from its generation and propagation to ultimate decay, but this is not essential to most experimental settings. Thus, a spike train can be described simply as a list of times, each indicating the time at which a given neuron emitted a spike. In statistics, this type of data is known as a point process: a spike can be considered a point in time.

Because a point is – at least in principle – infinitesimally small, no two spikes technically occur at exactly the same time. Yet, as the figure below shows, multiple presentations of the same stimulus (panel A) can evoke spike trains that are clearly similar across trials (panel B):

(from de Ruyter van Steveninck et al. 1997)

We would like a way to describe the observation that this neuron's response is similar across trials, even though the exact spike times differ between trials.

Spike counts (binned firing rate)

A simple idea is to count how many spikes occur in a given time interval (time bin). Then, the precise spike time is no longer relevant. Let's first plot a short snippet of an example spike train:

%% loading
cd 'C:\Users\mvdm\Dropbox\teaching\BIOL6xx\R042-2013-08-18'
fc = FindFiles('*.t');
S = LoadSpikes(fc);
%% plot
iC = 47;
t = [5801 5801.7];
spk_t = Data(Restrict(S{iC},t(1),t(2))); % get spike times
line([spk_t spk_t],[-1 -0.5],'Color',[0 0 0]); % note, plots all spikes in one command

To bin these spikes, we can use MATLAB's histc() function. This returns a histogram of spike counts, given a vector of bin edges:

binsize = 0.1;
tbin_edges = t(1):binsize:t(2); % vector of time bin edges (for histogram)
tbin_centers = tbin_edges(1:end-1)+binsize/2; % vector of time bin centers (for plotting)
spk_count = histc(spk_t,tbin_edges); % get spike counts for each bin
spk_count = spk_count(1:end-1); % ignore spikes falling exactly on edge of last bin.

When using histc, it is key to remember that element n of the output corresponds to the count in the interval [edge(n) edge(n+1)>. In other words, a spike occurring exactly at the time of the second edge will be assigned to output bin 2, not bin 1. As a result, histc() needs a special last bin for those spikes occurring exactly at the time of the last edge. The above code ignores this last bin, so that we have a number of output spike counts that is equal to the number of bin centers, which in turn is equal to the number of bin edges minus 1.

We can now plot the spike counts as follows:

hold on;
h = bar(tbin_centers,spk_count);
set(h,'BarWidth',1,'EdgeColor','none','FaceColor',[1 0 0]); % reformat bar appearance
yt = get(gca,'YTick'); ytl = get(gca,'YTickLabel');
set(gca,'YLim',[-1.5 10],'YTick',yt(2:end),'YTickLabel',ytl(2:end,:)); xlabel('time (s)');

This gives

☛ Verify by visual inspection that indeed the binned spike counts, shown as red bars, correspond to the number of spikes in each 100ms bin.

However, it should be clear that this count is a poor description of the underlying spike train. For instance, at time 5801.25 the spike count is zero, yet this time does not seem much different from, say, 5801.4. This occurs because we have sharp bin edges that are placed arbitrarily.

As an aside, note that histc() provides a raw spike count, which is not normalized by the bin size. Thus, to convert the counts into a firing rate estimate in spikes per second (Hz), we need to divide by the bin size.

Perhaps we can mitigate the problem of arbitrary edge placement by making the bin size smaller.

☛ Reduce the bin size to 10 ms and replot.

We now have a finer estimate that tracks the spike bursts and gaps more accurately, but it is still dissatisfying to see that some bins have a count of 2 and others of 1, implying that the firing rate for some bins was twice as large as those nearby, without clear justification. Again, how the bin edges happen to align with the spikes is the problem.

☛ What happens if you reduce the bin size to 1 ms?

At some point, if you make the spike count bins small enough, there will only ever be zero or one spikes per bin. This defeats the purpose of binning in the first place: now again each spike train will be unique. What we require is a measure that avoids the arbitrariness of bin edges.

A brief diversion: convolution

We saw in the previous section that binned spike counts run into problems both when bins are large (location of edges affect the count) and when bins are small (only one spike per bin at most, so no averaging). A way around this is to use small time bins to minimize the edge effects, combined with a certain amount of filtering (smoothing) to enable contributions from multiple spikes to a single bin.

As an example, imagine that each spike in the above spike train was replaced by a Gaussian, like this:

Then, we end up with a smoothly varying estimate of firing rate. Formally, this operation of replacing each spike (considered a Dirac delta function) with another function (in this example, a Gaussian) is known as convolution. Convolving a spike train is related to the idea of impulse response to a filter, which we encountered in an earlier module. Each spike is an impulse (input is zero everywhere except at the time of the spike) and the output is a filtered (convoluted) signal which is a smoothed version of the input.

Spike density function (SDF, binless firing rate)

To convolve our example spike train to obtain a smooth estimate of firing rate, we can use MATLAB's conv2() function:

binsize = 0.001; % select a small bin size for good time resolution
tbin_edges = t(1):binsize:t(2);
tbin_centers = tbin_edges(1:end-1)+binsize/2;
spk_count = histc(spk_t,tbin_edges);
spk_count = spk_count(1:end-1);
rw_sdf = conv2(spk_count,rectwin(50),'same'); % convolve with rectangular window
gau_sdf = conv2(spk_count,gausswin(50),'same'); % convolve with gaussian window

(The same option specified in the call to conv2() ensures that the output is the same size as the input.)

When plotted on top of the earlier binned spike train, you should get

Note how convolving with a rectangular function (in blue) results in an estimate that is jagged, but avoids the arbitrary bin edges imposed by the original spike count. For this estimate, each spike is effectively replaced with a rectangle centered at the spike. Convolving with a Gaussian (in green) avoids the jagged edges to make a smooth estimate.

Convolved spike trains are often referred to as spike density functions (SDFs) because they are effectively a continuous estimate of firing rate, i.e. describing the density of spikes over time. Of course, the SDFs we computed here are still discretized, but the resolution could be set to something arbitrarily fine by adjusting the bin size.

One final adjustment we need to make is to ensure we actually end up with a firing rate estimate in spikes per second (Hz). Because it is helpful to be able to state the convolution used in a succinct manner, it is best to use the gausskernel() function, which allows us to specify not only the width of the convolution window, but also the standard deviation of the Gaussian:

binsize = 0.001; % in seconds, so everything else should be seconds too
gauss_window = 1./binsize; % 1 second window
gauss_SD = 0.02./binsize; % 0.02 seconds (20ms) SD
gk = gausskernel(gauss_window,gauss_SD); gk = gk./binsize; % normalize by binsize
gau_sdf = conv2(spk_count,gk,'same'); % convolve with gaussian window

Note that we need to normalize by the binsize (1ms) because our firing rate estimate should be independent of the bin size chosen.

We thus obtain a SDF firing rate estimate with the correct units:

Of course, there is no single correct choice of the SD of the Gaussian kernel used for the convolution. This should be chosen by comparison with the raw spike train: in this case, we don't want to choose a Gaussian that is so wide that the modulation in firing rate is no longer visible.

It is also worth noting that convolution with a Gaussian to obtain a firing rate estimate has limitations. For instance, convolving spikes evoked by stimulus presentation might result in an increase in the firing rate estimate before the stimulus is presented. This can be undesirable and should not be interpreted as evidence of “stimulus anticipation” when in fact it simply results from the convolution! To minimize these sorts of issues, a variety of adaptive firing rate estimates have been developed, such as Bayesian Adaptive Regression Splines (BARS).

If you compute SDFs frequently, it makes sense to put the above into a function, which takes a ts object as input (along with binsize, windowsize and SD parameters), and returns a tsd firing rate estimate.

Peri-event or -stimulus time histograms

PETHs or PSTHs show the average firing rate around a time of interest, usually a task event such as the presentation of a stimulus. Once you have a firing rate estimate for the entire session, it is relatively straightforward to average this estimate over a fixed window around each event.

Interspike intervals (ISIs)

A different aspect of single neuron spike trains is the regularity of their spike timing. In principle, an average firing rate of say, 10 Hz could come about in a very regular manner (one spike exactly every 0.1s) or a more irregular manner (say, between 0.05 and 0.15s, with an average of 0.1). This idea is illustrated in the following figure, from a famous paper:

Panels A-C show simulated voltage traces (the membrane potential of a single model neuron), with the trace in panel A displaying very irregular firing, C very regular firing, and B something in between. How can we capture such differences?

ISI histograms and the coefficient of variation

A key difference between panels A-C above is the distribution of interspike intervals (ISIs). For a very regular spike train, there is a relatively stereotyped time that elapses before the next spike; for an irregular spike train, the interspike intervals are much more variable. Let's plot the ISI histogram for our own example neuron:

iC = 47;
spk_t = Data(S{iC}); % spike times
isi = diff(spk_t); % interspike intervals
dt = 0.001; % in s, because spike times are in s
isi_edges = 0:dt:0.25; % bin edges for ISI histogram
isi_centers = isi_edges(1:end-1)+dt/2; % for plotting
isih = histc(isi,isi_edges);
bar(isi_centers,isih(1:end-1)); % remember to ignore the last bin of histc() output
set(gca,'FontSize',20,'XLim',[0 0.25]); xlabel('ISI (s)'); ylabel('count'); grid on;

You should get

This histogram, which simply counts how many interspike intervals there are in this spike train for every interval from 0-250ms, shows that the most common interspike intervals are around 5ms. There is a hint of an increase in interspike intervals between 100-150ms, in the theta range typical for a hippocampal “place cell”. This ISI distribution can also be seen in the spike train snippet we used for computing firing rates, above: the short ISIs occur in bursts, and then there is a longer interval (in the theta range) until the next burst.

Clearly, this cell does not emit spikes at random intervals; instead, there is structure in the spike train. A further indication of this is the ISI return plot, a scatterplot of each interspike interval against the next.

☛ Using the isi variable created above, plot each ISI against its successor. Set the axes to run from 0 to 0.25 seconds.

You should get something like this:

This plot shows that a short ISI (in the 5ms range) tends to be followed either by another short ISI, or by a long one in the theta range. It however seems rare for a theta-range ISI to be followed by another theta-range ISI, again in accordance with our visual impression of the spike train.

Interlude: Poisson point process

A useful tool for many spike train analysis questions is to make a comparison with a random, synthetic spike train. Such comparisons can bring into focus salient aspects of real spike trains, and can be used for tests of significance.

A simple way of generating a random spike train is to essentially flip a coin at each time step, scoring a spike for heads, and and no spike for tails:

dt = 0.001;
t = [0 10]; % time interval (length) of spike train to generate
tvec = t(1):dt:t(2);
pspike = 0.5; % probability of generating a spike in bin
rng default; % reset random number generator to reproducible state
spk_poiss = rand(size(tvec)); % random numbers between 0 and 1
spk_poiss_idx = find(spk_poiss < pspike); % index of bins with spike
spk_poiss_t = tvec(spk_poiss_idx)'; % use idxs to get corresponding spike time
line([spk_poiss_t spk_poiss_t],[-1 -0.5],'Color',[0 0 0]); % note, plots all spikes in one command
axis([0 0.1 -1.5 5]); set(gca,'YTick','');

This gives:

This spike train looks different from our real one.

☛ Plot its ISI histogram. Some differences are apparent; for instance, this synthetic spike train clearly has much shorter ISIs overall.

This difference occurs because the mean firing rates between the two spike trains are not equal. Our real spike train has a mean firing rate of about 0.47 spikes/s (as 1./mean(diff(spk_t)) will verify). In contrast, our synthetic spike train has a theoretical mean firing rate of 500 spikes/s (50% probability of a spike in each 1ms bin).

☛ What should the probability per 1ms bin of a spike be, in order for our synthetic spike train to have a mean firing rate of 0.47 spikes/s? Recompute the ISI histogram for a spike train generated with this probability instead of 0.5. You will need to increase the length of the spike train to 4500s in order to get a number of spikes which is similar to that in the real spike train.

You should get something like:

The number of “short” ISI counts (i.e. the ones shown in this plot) is way smaller in our synthetic spike train than in the real one. So, even though the mean firing rates are the same, our real spike train has relatively many short ISIs, which must be balanced with longer ones to end up with the same mean.

Generating artificial spike trains in this way, by flipping a coin at each time step, is known as a Poisson point process (Poisson process for short): it generates interspike intervals drawn from the Poisson distribution. An important property of this distribution is that its values are never negative. This is useful when dealing with counts and interspike intervals, which cannot be negative. For large means, the Poisson distribution will approach a Gaussian distribution, but a Gaussian with a mean close to zero is likely to generate negative values, which are nonsensical when dealing with ISIs or spike counts.

A further property of the Poisson process that each ISI is independent of the last.

☛ Generate the ISI return plot for your Poisson spike train. If you only see a few spikes, increase the axes limits.

You should see a total absence of structure (correlations) in the plot, unlike the equivalent for the real spike train, indicating that a given ISI does not tell you anything about what the next one might be.

A final difference with our actual spike train is that pure Poisson spike trains have no refractory period. Verify that the ISIh of the real spike train does!

Our synthetic Poisson spike train had a fixed probability of emitting a spike in each time bin. This results in a homogenous Poisson process, meaning that theoretically, the firing rate does not vary over time. Of course, there will be random fluctuations, but the underlying probability of a spike remains constant. In contrast, an inhomogenous Poisson process involves changes over time of the underlying spike probability.

☛ In the above code for generating a spike train, pspike was set to a fixed probability of 0.5. Change pspike to be a vector of the same size as tvec describing a 1 Hz sinusoid which oscillates between 0 and 0.5. Use this changing spike probability to generate your artificial (inhomogenous Poisson) spike train. Verify by visual inspection that indeed the firing rate of this spike train is varying over time.

Spike autocorrelation function (acf)

A different way of characterizing the properties of a spike train is to compute its autocorrelation function, which describes the probability of finding a spike at time $t+t'$ given a spike at time $t$, for some range of lags $t'$. This is analogous to the autocorrelation of gamma-band power in a ventral striatal LFP that we computed in the previous module.

MATLAB doesn't provide a built-in function that can take in a spike train and compute its autocorrelation, so we have to make our own. Here is a simple implementation:

function [ac,xbin] = acf(spike_times,binsize,max_t)
% function [ac,xbin] = acf(spike_times,binsize,max_t)
% estimate autocorrelation of input spike train
% spike_times: [1 x 1 ts]
% binsize: acf bin size in s
% max_t: length of acf in s
% ac: autocorrelation estimate (normalized so that acf for the zero bin is
% 0)
% xbin: bin centers (in s)
% MvdM 2013
spk_t = Data(spike_times);
xbin_centers = -max_t-binsize:binsize:max_t+binsize; % first and last bins are to be deleted later
ac = zeros(size(xbin_centers));
for iSpk = 1:length(spk_t)
   relative_spk_t = spk_t - spk_t(iSpk);
   ac = ac + hist(relative_spk_t,xbin_centers); % note that hist() puts all spikes outside the bin centers in the first and last bins! delete later.
xbin = xbin_centers(2:end-1); % remove unwanted bins
ac = ac(2:end-1);
ac = ac./max(ac); % normalize

☛ Using this autocorrelation function and the code for generating Poisson spike trains above, compute and plot the autocorrelation of a Poisson spike train with a mean firing rate of 0.47 Hz (exactly as done above), a 10ms time bin, from -1 to 1 seconds. Note that the acf() function requires a ts object as input, so you will need to convert the spike times to a ts object using ts(spk_poiss_t).

You should get something like:

The autocorrelation of a Poisson spike train appears to be essentially flat over different time lags.

☛ Why? What differences to you notice with the acorr of the actual spike train ([acorr,xbin] = acf(S{47},0.01,1);)?

Spike cross-correlation function (ccf)

The intuition behind the cross-correlation between two spike trains is the following: given that neuron 1 emits a spike at time $t$, what is the probability of seeing a spike from neuron 2 at different lags $t+t'$? For instance, if neuron 2 reliably fires 10ms after neuron 1, the cross-correlation between 1 and 2 would have a peak at +10ms. The autocorrelation is a special case of this with the two spike trains being the same.

Computing the cross-correlation is a simple modification of our autocorrrelation function:

function [cc,xbin] = ccf(spike_times1,spike_times2,binsize,max_t)
% function [cc,xbin] = ccf(spike_times1,spike_times2,binsize,max_t)
% estimate cross-correlation function of input spike train
% spike_times1: [1 x 1 ts]
% spike_times2: [1 x 1 ts]
% binsize: ccf bin size in s
% max_t: length of acf in s
% cc: cross-correlation estimate (relative to spike_times1)
% xbin: bin centers (in s)
% MvdM 2013
spk_t1 = Data(spike_times1);
spk_t2 = Data(spike_times2);
xbin_centers = -max_t-binsize:binsize:max_t+binsize; % first and last bins are to be deleted later
cc = zeros(size(xbin_centers));
for iSpk = 1:length(spk_t1)
   relative_spk_t = spk_t2 - spk_t1(iSpk);
   cc = cc + hist(relative_spk_t,xbin_centers); % note that histc puts all spikes outside the bin centers in the first and last bins! delete later.
xbin = xbin_centers(2:end-1); % remove unwanted bins
cc = cc(2:end-1);
cc = cc./(length(spk_t1)); % normalize by number of spikes of first input

The only potentially tricky part of cross-correlations is how to do the normalization (there are some different choices that determine if the output is a conditional probability, a correlation, or something else), but we will not worry about this for now. All these normalizations produce outputs that are qualitatively similar.

☛ What do you expect the cross-correlation function of two Poisson spike trains to look like? Why? Verify your intuition.

Let's compute the cross-correlation between two hippocampal place cells, simultaneously recorded from a rat running a T-maze:

cell1_id = 5; cell2_id = 42;
s1 = Restrict(S{cell1_id},3200,5650); % restrict to on-track times only
s2 = Restrict(S{cell2_id},3200,5650);
[xcorr,xbin] = ccf(s1,s2,0.01,1);
set(gca,'FontSize',20); xlabel('lag (s)'); ylabel('xcorr');

You should get:

This is an example of one of my favorite plots in neuroscience. It holds the key to revealing a fundamental and unique property of the rodent hippocampus, thought to reflect a specialization for the rapid encoding of episodic memories. The plot shows (something akin to) the probability of cell 42 spiking at various time lags (between -1 and 1 second) given that cell 5 spikes at time 0. This ccf is clearly asymmetric, in that cell 42 is more likely to spike after cell 5 than before, as evident from the crosscorrelation values on the right side of the plot (5–>42) being larger overall than those on the left (42–>5).

There is something else going on, too: on top of this slow component of the ccf, there is a faster component with sharp peaks repeating at theta frequency (about 8-9 per second). The peak closest to zero is in fact slightly offset, occurring at approximately 30-40ms. Thus, cell 42 is likely to spike just a few tens of ms after cell 5; in contrast, the reverse, where cell 5 spikes a few tens of ms after cell 42, almost never happens.

As it turns out, this highly precise ccf shape results from the fact that each theta cycle in the hippocampal LFP is in fact associated with the sequential activation of a number of place cells, tracing out a trajectory on the maze, as illustrated in this figure from Skaggs et al. 1996:

A different way of showing this idea is as follows (from Malhotra et al. (2012):

This figure illustrates the ccf (labeled 'CCG' here, for cross-correlogram, which is the same thing) that would be expected from generating inhomogenous Poisson spike trains based on a firing rate profile given by two place fields (the blue and red lines in the top panel). In this case, with the animal running from left to right, the blue spikes would tend to come after the red, because the red field is entered before the blue field. This is reflected in the asymmetric shape of the ccf. However, this lacks the sharp theta peaks apparent in the ccf we just computed from our two real place cells. The figure here shows that if we enable an effect called “phase precession” (see here for more optional info on this) then the shape of the ccf does reproduce what is found experimentally. You can also see that accordingly, the blue spikes tend to come just after the red spikes within each theta cycle, forming sequences of place cells.

An important question is to establish if such sequences could come about by just having each place cell independently emit spikes in a manner described by some average spike density function, or if the sequences that are actually seen require some sort of coordination between place cells.

To test this, we can first convert our spike trains to a SDF over time, and then use that SDF to generate inhomogenous Poisson spike trains. Since these two spike trains are generated by independent coin flips, i.e. whether or not there was a spike for neuron 1 cannot directly affect the probability of a spike in neuron 2, we can determine if this is sufficient to reproduce the ccf seen in the real data (in which there may be some dependency between neuron 1 and neuron 2).


☛ Compute the ccf of two spike trains generated from the SDF of our two hippocampal place cells (5 and 42) and compare it to their actual ccf. Use a 50 ms standard deviation for the Gaussian, and 1ms bin size. The sequence of steps for this would look something like this:

% for each of the two neurons, restrict the data to [3200 5650] (the time interval when the rat was running on the track)
% compute the spike density function for each, making sure that your tvec runs from 3200 to 5650 also, and that you have a 50ms SD for the Gaussian convolution kernel
% to use these SDFs to generate Poisson spike trains, convert the firing rates given by the SDF to a probability of emitting a spike in a given bin. (As you did above for a 0.47 Hz constant firing rate.)
% generate Poisson spike trains, making sure to use the same tvec
% convert Poisson spike trains to ts objects and compute the ccf

(This comparison of independently generated spikes with the actual data is the essence of this excellent paper.)


Erwin, 2019/11/29 12:07

I love browsing your internet site. thnx! Generic viagra 100mg canadian women generic viagra sildenafil citrate reviews vs viagra generic cheap generic viagra onlin pharmacy in south carolina 2016 [url=]sildenafil generic viagra reviews[/url] lady era female viagra reviews youtube

Thomasblapy, 2019/12/14 12:08

canadian pharmacy viagra <a href=“ ”>walmart online pharmacy</a> walmart pharmacy online

the pharmacy, 2019/12/15 02:25

the canadian pharmacy modafinil online pharmacy

GeorgeSpifT, 2019/12/15 07:10

canada drugs <a href=“ ”>legitimate canadian pharmacy</a> vipps canadian pharmacy

Kevindiady, 2019/12/15 10:18

canadian pharmacy discount code <a href=“ ”>24 hour pharmacy near me</a> legit canadian pharmacy <a href=>24 hr pharmacy</a> pharmacy discount card

GeorgeSpifT, 2019/12/15 12:27

california pharmacy <a href=“ ”>canadian pharmacy cialis reviews</a> cheapest pharmacy

Kevindiady, 2019/12/15 16:20

cialis online pharmacy <a href=“ ”>indian pharmacy online</a> canadian pharmacy world coupons <a href=>precription drugs from canada</a> canadian pharmacy cialis

GeorgeSpifT, 2019/12/15 17:44

legit online pharmacy <a href=“ ”>modafinil online pharmacy</a> canadian pharmacy near me

Kevindiady, 2019/12/15 22:06

the canadian pharmacy <a href=“ ”>viagra canadian pharmacy vipps approved</a> canadian pharmacy levitra value pack <a href=>online canadian pharmacy review</a> legit canadian pharmacy

GeorgeSpifT, 2019/12/15 23:03

canadian drug stores <a href=“ ”>pharmacy discount card</a> legit canadian pharmacy online

Kevindiady, 2019/12/16 04:07

best canadian online pharmacy reviews <a href=“ ”>vipps canadian pharmacy</a> medical pharmacy <a href=>canadian pharmacy online cialis</a> online canadian pharmacy

GeorgeSpifT, 2019/12/16 04:38

online pharmacy <a href=“ ”>canadian mail order pharmacy</a> pharmacy online

AndrewFak, 2019/12/16 15:41

canadian viagra <a href=“ ”>viagra canada</a> cheap viagra <a href=>viagra without a doctor prescription walmart</a> cialis vs viagra

GlennfiT, 2019/12/16 15:42

generic viagra 100mg

GlennfiT, 2019/12/16 20:58

viagra substitute over counter

AndrewFak, 2019/12/17 03:18

sildenafil citrate generic viagra 100mg <a href=“ ”>buy viagra online</a> female viagra <a href=>viagra coupons 75% off</a> viagra cost per pill

GlennfiT, 2019/12/17 06:23

generic viagra prices

GlennfiT, 2019/12/17 12:18

buy viagra online

GlennfiT, 2019/12/17 18:08

viagra samples

GlennfiT, 2019/12/18 00:04

generic viagra available

AndrewFak, 2019/12/18 02:33

viagra without doctor prescription <a href=“ ”>sildenafil vs viagra</a> viagra price <a href=>generic viagra</a> what does viagra do

GlennfiT, 2019/12/18 05:50

how long does viagra last

GlennfiT, 2019/12/18 17:29

how long does viagra last

AndrewFak, 2019/12/18 19:42

cost of viagra 100mg <a href=“ ”>herbal viagra</a> natural viagra <a href=>viagra cost per pill</a> buy viagra

Steveninsiz, 2019/12/19 08:46

viagra pills <a href=“ ”>how to take viagra for maximum effect</a> when will viagra become generic

Fletcher, 2020/06/27 08:28

male enhancements viagra online viagra professional

viagra online, 2020/06/27 08:28

male enhancements viagra online viagra professional

天天领现金, 2020/09/04 05:36


天天领现金, 2020/09/06 09:06


GScraper, 2020/09/08 07:43


GScraper, 2020/09/10 05:38


5 Inch Usb Fan, 2020/09/10 23:54

5 Inch Usb Fan

1x22 Red Dot Sight, 2020/09/11 00:12

1×22 Red Dot Sight

0.5 Mm Plastic Sheet, 2020/09/11 01:50

0.5 Mm Plastic Sheet

Electric Pressure Sustaining Control, 2020/09/11 03:11

Electric Pressure Sustaining Control

Best High Bay Lights, 2020/09/11 03:49

Best High Bay Lights

Gabion Mesh Machine, 2020/09/11 07:12

Gabion Mesh Machine

Boiler Equipment, 2020/09/11 07:33

Boiler Equipment

Auxiliary equipment, 2020/09/11 07:49

Auxiliary equipment

Brake drum, 2020/09/11 09:08

Brake drum

Checkweigher System, 2020/09/11 09:31

Checkweigher System

N-Octyl-2-Pyrrolidone, 2020/09/11 11:04


Sporting Type Scope, 2020/09/11 12:13

Sporting Type Scope

Galaxy Black Marble, 2020/09/11 12:48

Galaxy Black Marble

17mm Pvc Foam Board, 2020/09/11 13:20

17mm Pvc Foam Board

0.2mm Thick Pvc Rigid Plastic Sheet, 2020/09/11 14:42

0.2mm Thick Pvc Rigid Plastic Sheet

0445120007, 2020/09/11 15:10


Blue Marble Mosaic, 2020/09/11 15:35

Blue Marble Mosaic

Aluminium Lift Sliding Glass Doors, 2020/09/11 16:02

Aluminium Lift Sliding Glass Doors

Statures, 2020/09/11 16:56


Best Led Flat Panel Lights, 2020/09/11 17:25

Best Led Flat Panel Lights

Check Vavle, 2020/09/11 17:53

Check Vavle

15ml airless pump bottles, 2020/09/11 19:29

15ml airless pump bottles

Chain Link Fencing Die Machine, 2020/09/11 19:56

Chain Link Fencing Die Machine

Bamboo Sheet Sets, 2020/09/11 20:04

Bamboo Sheet Sets

Pillow Case, 2020/09/11 20:45

Pillow Case

1000l Beer Brewing House, 2020/09/11 20:58

1000l Beer Brewing House

Barbed Wire Making Machine, 2020/09/11 21:19

Barbed Wire Making Machine

FTTx Equipment, 2020/09/11 23:10

FTTx Equipment

Airsoft Scope Mount, 2020/09/12 00:48

Airsoft Scope Mount

25G SFP28 to SFP28, 2020/09/12 01:38

25G SFP28 to SFP28

1/25 4mm Alum Rings, 2020/09/12 03:55

1/25 4mm Alum Rings

Granite Hollow Columns, 2020/09/12 05:19

Granite Hollow Columns

4x8 Foam Sheets, 2020/09/12 05:48

4×8 Foam Sheets

Discontinued Marble Tile, 2020/09/12 07:31

Discontinued Marble Tile

Chevron Pattern Mosaic Tile, 2020/09/12 13:32

Chevron Pattern Mosaic Tile

0 928 400 627, 2020/09/12 15:30

0 928 400 627

Brown Marble Mosaic Tile, 2020/09/12 15:57

Brown Marble Mosaic Tile

Cappuccino Marble Mosaic, 2020/09/12 21:04

Cappuccino Marble Mosaic

Anti Static Plastic Film, 2020/09/12 22:26

Anti Static Plastic Film

China Gazebo Sculpture, 2020/09/12 23:47

China Gazebo Sculpture

Despirt Mosaic & Marble, 2020/09/13 05:10

Despirt Mosaic & Marble

American Style, 2020/09/13 07:07

American Style

Glass Pebble Mosaic Tile, 2020/09/13 07:28

Glass Pebble Mosaic Tile

Body Sensor Heater, 2020/09/13 07:49

Body Sensor Heater

1mw Laser Sight, 2020/09/13 09:11

1mw Laser Sight

100G QSFP28 to 4x SFP28, 2020/09/13 11:11

100G QSFP28 to 4x SFP28

Adjustable Scope Base, 2020/09/13 11:48

Adjustable Scope Base

Chicken Packaging Machine, 2020/09/13 13:21

Chicken Packaging Machine

Electric Section, 2020/09/13 15:18

Electric Section

100 Gallon Beer Ferment, 2020/09/13 17:07

100 Gallon Beer Ferment

Welded Wire Mesh Machine, 2020/09/13 17:33

Welded Wire Mesh Machine

100ml lotion bottles, 2020/09/13 19:18

100ml lotion bottles

100w Ufo High Bay, 2020/09/13 19:36

100w Ufo High Bay

Aluminium Checkered Sheet 3003, 2020/09/13 21:36

Aluminium Checkered Sheet 3003

Volvo Injector, 2020/09/13 23:02

Volvo Injector

biomass boiler water tank, 2020/09/14 03:32

biomass boiler water tank

Artemide Linear, 2020/09/14 08:48

Artemide Linear

Armored Optical Cable, 2020/09/14 13:15

Armored Optical Cable

CV joint kit, 2020/09/14 13:52

CV joint kit

boiler accessory slag remover, 2020/09/14 15:24

boiler accessory slag remover

1000l Brewery, 2020/09/14 15:53

1000l Brewery

3 Way Cock Valve- Brass, 2020/09/14 17:49

3 Way Cock Valve- Brass

Coffee Tables, 2020/09/14 19:13

Coffee Tables

Printing base film, 2020/09/14 19:52

Printing base film

Attacco su Titano Fairy Tail Tokyo Ghoul Canvas Horizontal Square Borsa A Tracolla. Borsa A Tracolla. Fantasia A Fumetti. Usata per Cosplay. Fumetti, 2020/09/15 01:58

Bolso de playa de verano para mujer Bolso de paja (Beige) Damen Sandalen Plateau Keilabsatz Offener Zeh Espadrilles Frauen Kn?chel Schnalle Strandschuhe.06.39 Sud Trading Company. Etui permis de Conduire Fée pachier Aluminium Frame Insulated Windows Attacco su Titano Fairy Tail Tokyo Ghoul Canvas Horizontal Square Borsa A Tracolla. Borsa A Tracolla. Fantasia A Fumetti. Usata per Cosplay. Fumetti

Chaussures de sécurité basses Dickies Norden II S3 SRC, 2020/09/15 07:12

Automatic Chain Link Fencing Machine City Flarelo T. Mocasines para Ni?os Damen Beacon V3 Fresh Foam Laufschuh Signature F16109 - Portafoglio in PVC con doppia cerniera. colore: Cachi Chaussures de sécurité basses Dickies Norden II S3 SRC

Frauen über den Knie Stiefel Mode Winter warm halten Schneeschuhe Mode M?dchen Bequeme Schuhe, 2020/09/15 11:13

Star Runner 2. Zapatillas de Atletismo para Ni?os Fiberhome ONU GirlsDrew Barrymore Isabella Sandal Pantofole infradito da donna con zeppa in PU (poliuretano) comfort estivo nero/argento/rosa. argento.US6.5-7 / EU37 / UK4.5-5 / CN37 Frauen über den Knie Stiefel Mode Winter warm halten Schneeschuhe Mode M?dchen Bequeme Schuhe

Monedero para mujer. cierre de clip. 100% piel de anguila. Carrot (Naranja) - 100004-642, 2020/09/15 11:55

Femmes Bout Rond Escarpins Plate-Forme Chaussures habillées Sexy Mode fête de Mariage Talons Aiguilles en Cuir Verni Talon Mince Talons Hauts biomass boiler induced draft fan Damen Charm Sunflower Handbag Shz 1 Henkeltasche. Blau (Light Blue). 16x20x28 cm 7160-665. Portafogli Donna CUOIO talla unica Monedero para mujer. cierre de clip. 100% piel de anguila. Carrot (Naranja) - 100004-642

0.3mm Rigid Transparent Pvc Rolls, 2020/09/15 13:33

1-10 8004. Sandalias con Plataforma para Mujer Neivobos Frauen ethnische Baumwolle und Leinen Eimer Tasche gewebt Umh?ngetasche Retro Quaste Umh?ngetasche (Color : Green. Size : 23 * 18 * 27cm) Fox Geldb?rse The Corner Chuck Taylor all Star Galaxy Shimmer High Top 0.3mm Rigid Transparent Pvc Rolls

Northeast Agricultural University, 2020/09/15 15:05

Sacs à main pour femmes chien dr?le lecture journal fourre-tout en cuir PU Moschino kids Borse Neonato Rosso Mqx036 lda16 poppy red Unisex Zapatos De Vadeo Botas De Monta?a Enviar Calcetines Unisex Kinder Crocsfl Paw Patrol Band K Clogs Northeast Agricultural University

Edge Lux 3 W. Chaussure de Course Femme, 2020/09/15 15:35

Women Small Cell Phone Purse Crossbody.Colorful Antlers Boho Deer Retro Artsy Winter Rain Pattern Rainbow Heart Animal Theme Notting Hill Cartera Minimalista con Cremallera y protección RFID para Tarjetas Monedas y Efectivo Color Negro Winterstiefel Damen Gefüttert Wasserdicht Winterschuhe Warm Boots Winter Kurzschaft Stiefel für Damen Flach Stiefeletten rutschfeste Wanderschuhe Schwarz blau rot Building Concrete Mesh Welding Machine Edge Lux 3 W. Chaussure de Course Femme

Summer II 8 Velour. Scarpe con Lacci Uomo. Black. 42.5 EU, 2020/09/15 17:36

Stellar Tac Botas militares y tácticas impermeables para hombre Chaussons de Bain Douche Sandales Corgi Chaussures toboggans de Plage Sandale antidérapante Zipperb?rse Mercy. 15 cm. 11985 Accurate Laser Bore Sighter Summer II 8 Velour. Scarpe con Lacci Uomo. Black. 42.5 EU

Air First Class. Scarpe da Fitness Uomo, 2020/09/15 19:33

Damen Babson Stiefelette Flowmeter KUCCI Zapatos de Vestir Tacón Alto Clásico Moda Punta Puntera para Mujer Women Small Cell Phone Purse Crossbody.Creative Composition With Dots And Stripe Monochrome Design Brush Mark Effect Air First Class. Scarpe da Fitness Uomo

Sacs à Bandoulière Pour Femmes Grand Sac Hobo Bucket Purse Toile Sac à Bandoulière Sac à Bandoulière.Khaki, 2020/09/15 19:52

Damen Cai Umh?ngetasche. 28x22x14 cm 10 Bbl Brite Tank For Beer Remaches de cuero de las mujeres de la parte superior de la manija de las bolsas cruzadas del cuerpo de la moda mochilas para las compras trabajo campus negro calabaza boceto Wguili Uomo Sandali Mens Sport Outdoor Sandali Chiuso Le Dita dei Piedi Scarpe di Cuoio Estate Pescatore Beach per Camminare (Color : Black. Size : 42) Sacs à Bandoulière Pour Femmes Grand Sac Hobo Bucket Purse Toile Sac à Bandoulière Sac à Bandoulière.Khaki

CafePress Doctor Nurse Sac fourre-tout en forme de c?ur. Toile. kaki. Taille M, 2020/09/16 05:16

Derivatives?Resource Sandalias Alpargatas Plataforma Tacón De Cu?a Cordón Arriba para Mujeres Donna Stivaletti Bassi Stivaletti Stringati 144216 Nero 39 Flandell Alwayswin Herren Winterstiefel Plus Samt Warme Kurze Stiefel Outdoor rutschfeste SchneeschuheTrekking-Wanderstiefel Leder wasserdichte Stiefel Stiefeletten Schnürstiefel Turnschuhe Booties CafePress Doctor Nurse Sac fourre-tout en forme de c?ur. Toile. kaki. Taille M

HUOQILIN Taschen Handtaschen Mode Wilde Breitband Eimer Tasche Umh?ngetasche Messenger Bag (Color : Red), 2020/09/16 07:26

East China University Of Science and Technology Uomo Sandali Pelle Sandalo Uomini Morbida Outdoor Closed Toe Shoes Casual Sandali Sport Spiaggia Ballerine Mano Stitching Slip On Wide-Fit Shoes Pescatore di Lavoro Adatto ad Usura di Estate Secado Rápido Volibear CRCOG Hommes Cikou Peu Court de la Carte Porte-Monnaie Paragraphe rétro Portefeuille Mince Moins Cher CRCOG (Color : Light Brown. Size : 14X10CM) HUOQILIN Taschen Handtaschen Mode Wilde Breitband Eimer Tasche Umh?ngetasche Messenger Bag (Color : Red)

Scarpe Basse Sneakers Donna Rosa (DISRUPTOR-2-EMBROIDERY_5FM00605), 2020/09/16 07:47

Wmns Quest. Scarpe Running Donna Chain Link Jali Machine ARIAT Damen Heritage Iv French Paddock Boot Ooahh Sport. Pantuflas Unisex Adulto Scarpe Basse Sneakers Donna Rosa (DISRUPTOR-2-EMBROIDERY_5FM00605)

KROSER Laptop Damen Handtasche Shopper 15.6 Stilvolle Umh?ngetasche Wasserabweisende Gro?e Reise Einkauftasche mit RFID-Taschen für Arbeit/Business/Hochschule/Frauen-Tarnung Schwarz MEHRWEG, 2020/09/16 11:13

TIZORAX Boston Terrier - Borse per scarpe da viaggio. impermeabile. per uomo e donna Chengxin Botas de Combate de Media Pantorrilla Botas Altas de Hombres. ata for Arriba del Cuero auténtico Suela de Goma de Costura de Doble Correas Monk Shoes (Lana Forrada Opcional) Blue Onyx Pattern Rose Love Girl médite Peace Dove Vintage Pouch Girl -Lock Changement Sac à Main Portefeuilles Boucle en Cuir Porte-Monnaie clé Femme imprimé KROSER Laptop Damen Handtasche Shopper 15.6 Stilvolle Umh?ngetasche Wasserabweisende Gro?e Reise Einkauftasche mit RFID-Taschen für Arbeit/Business/Hochschule/Frauen-Tarnung Schwarz MEHRWEG

Reflector base film, 2020/09/16 11:58

Bobs - Ballet plano para mujer Champions Are The Breakfasts of Chuck Norris - Zaino con porta di ricarica USB Artnr. 7548 DUNLOP PUROFORT MULTI GRIP SAFETY Stiefel. PU. blau. mit Stahlkappe. erh?ltliche Gr??en: 36-49 Harvesthouse Beautiful Sea Under Blue Sky and White Clouds Custom?Wooden?Jigsaw?Puzzles?Kids?Adults?Toys?DIY?Wooden?Jigsaw?Puzzles?Family?Wall?Decoration Reflector base film

Perfume Spray Pen, 2020/09/16 13:35

Borse piccole a tracolla Borsa per cellulare Piccola stampa floreale con fessure per carte di credito TuToy Einstellbare Stiefelspanner Breite Schuhformer Kiefer Holzstiefel Baum Stretch Für M?nner Frauen Eu35-46 - # 2 Sac en Toile Sac à provisions réutilisable Simple Sac Quotidien décontracté Sac à Main Sac à bandoulière Remaches de cuero de las mujeres de la parte superior de la manija de las bolsas cruzadas del cuerpo de la moda mochilas para el trabajo de compras campus negro rico sopa de verduras Perfume Spray Pen

Ultrarange Exo Se EU 42, 2020/09/16 15:06

Products Uomo Donna Scarpe da Ginnastica Corsa Sportive Fitness Running Palestra Sneakers Basse Scarpe Comode per Camminare Jogging Nstqlzh Freizeit Driving Loafers for M?nner beil?ufige Flache Schuhe Penny-runde Zehe weicher Mikrofaser Leder Perforierte Beleg auf Lightweight` (Color : Brown-Breathable. Size : 37 EU) Grand sac de courses Tunis pour femme Ultrarange Exo Se EU 42

C.T. All Star Hi LTD Sneakers TESSUTO WHITE SMOKE GRIGIO 156885C, 2020/09/16 15:59

Electric Steam Boiler Gep?ck Handtaschen Twill Hoop Handtasche Stereo Paket Women Small Cell Phone Purse Crossbody.Pattern Of Chrysanthemum Flowers With Retro Inspirations Flourishing Nature Design Tongs ergonomiques Orteils Tongs.Sandales à Talons compensés et Pantoufles avec Strass.Chaussures de Plage-Neuf Rouge_40.Chaussures de Plage et de Piscine pour la Douche C.T. All Star Hi LTD Sneakers TESSUTO WHITE SMOKE GRIGIO 156885C

Official Star Wars Flap Purse, 2020/09/16 17:25

Borsa a tracolla in pelle rotonda riccia albero carino moda donna borsa a tracolla regolabile borsa superiore regolabile per donna ragazza HROIJSL Herren im britischen Stil mit Freizeitschuhen und Business-Schuhen Herrenmode L?ssig Feste Schnürschuhe Leder Hochzeitsschuhe M?nnliche Herrenschuh Osmond Schlichter Business-Halbschuh 0.2mm Packing Pvc Sheet Porte-Monnaie Nouveautés Trifold Women Wallet Titulaire de la Carte Porte-Monna