My Scipy code that generates data, plots it and performs the fft is listed below. The code generates two "spike trains" at intervals of 6 and 8 seconds. The following plot shows the signal in time, and its fft (sorry, the fft title got mangled).

The interesting thing to notice is that the frequencies you are looking for are there. Said frequencies are
1/8 Hz or omega = 1/8/500 = 0.00025 omega and
1/6 Hz or omega = 1/6/500 = 0.00033
Which is where the first two peaks are in the frequency domain. Further out are many more repeating peaks/aliases because these are essentially really extreme, intermittent square waves. You could try low-pass-filtering to isolate the fundamental frequencies and eliminate other peaks, but I think that's probably a lost cause at a sample rate of 500 Hz. Trying to get a 0.2 Hz cutoff filter @ 500 Hz seems a bit much. Just zooming in on the frequency region of interest might be good enough. The next step could be to understand exactly how these spike trains are going to cause repetition/aliasing in the frequency domain and look for a certain pattern.
Or try to either LPF and downsample a lot and retry the FFT route, calculate individual DFT bins using straight-up DFT math or use the Goertzel algorithm to calculate DFT bins for suspected frequency locations. Maybe try a Goertzel algorithm from 0.0 Hz to 0.2 Hz with 0.01 spacing - just as an example.
Good luck - it doesn't seem like a trivial problem.
from pylab import *
close('all') #close previous plots
fs=500 # generate signal
s1 = hstack(( 1,zeros(fs*6-1) )) # 1/6 Hz
s2 = hstack(( 1,zeros(fs*8-1) )) # 1/8 Hz
s1r = tile(s1,(12,))
s2r = tile(s2,(9,))
s = s1r + s2r
print "shape of s is %s"%(str(s.shape))
sfft=fft(s) # execute fft
sfreqs=fftshift(fftfreq(len(sfft)))
print "shape of sfft is %s"%(str(sfft.shape))
figure(); # generate plots
subplot(211); plot( 1./fs * arange(len(s)) , s ) #, 'o' )
title('signal, s');xlabel('sec');ylabel('mag')
subplot(212); plot( sfreqs , (abs(sfft)) ) #, 'o' )
title('fft(s)');xlabel('$\omega$');ylabel('linear mag')
axis([0.0,0.0005,0.0,25.])
show()
UPDATE: Thinking about it some more, the best solution probably doesn't involve signal processing. You may want to threshold the samples in to zeros and spikes. Then take the interval between the spike locations 0 and 1 and look for that interval through the rest of the spikes. If that interval is consistent with all of your spikes, then you've found that series of spikes and should "subtract out" that series from the spikes. Then keep iteratively looking at intervals of spikes. In a noise, imperfect system, a good implementation of this will probably be tricky, but it's probably a better solution to this specific problem.