% synapse_poisson.m generates Poisson trains of spikes with changes in rate
% then calculated synaptic changes generated by the spike train 
clear
dt = 0.0001;                            % time step
tmax=10;                                % maximum time

p0 = 0.5;                               % base vesicular release probability
tau_syn2 = 0.002;                       % determines synaptic rise time
tau_syn1 = 0.100;                       % synaptic decay time
tau_d = 0.5;                            % time for recovery from depression
tau_f = 0.5;                            % time for recovery from facilitation
f_f = 0.0;                              % facilitation factor: increases vesicle release probability
tau_syn = max(tau_syn1,tau_syn2);       % synaptic time constant used for calculation

t=dt:dt:tmax;                           % series of time points
Nt = tmax/dt;                           % number of time points

Nrates = 4;                             % number of different average rates

ratevalue(1)=25;                        % initial firing rate
ratevalue(2)=100;                       % second firing rate
ratevalue(3)=10;                        % third firing rate
ratevalue(4)=40;                        % fourth firing rate

rate = zeros(size(t));                  % array to contain rate as a function of time
for n = 1:Nrates                        
    for i = 1:Nt/Nrates                         
        rate(i+(n-1)*Nt/Nrates) = ratevalue(n); % puts appropriate rate at the right time
    end
end

% First calculate the expected response for a depressing synapse

Dmean = zeros(size(t));                 % initialize vector of calculated mean values for D
Synmean = zeros(size(t));               % initialize vector of calculated mean valued for synaptic s

Dmean(1) = 1.0;                         % start with no depression before spikes
Synmean(1) = 0.0;                       % start with no synaptic conductance open  
for i = 2:length(t)
    Dmeaninf = 1/(1+rate(i)*p0*tau_d);      % calculated mean depression rate for a Poisson train
    tau_deff = tau_d/(1+Dmean(i-1)*rate(i)*tau_d);    % time constant for depression depends on rate 
    Dmean(i) = Dmeaninf - (Dmeaninf-Dmean(i-1))*exp(-dt/tau_deff);  % update calculated value of D
    
    Synmeaninf = p0*Dmean(i)*rate(i)*tau_syn/(1+p0*Dmean(i)*rate(i)*tau_syn); % use this to estimate 
                                                                              % synaptic conductance
    Synmean(i) = Synmeaninf - (Synmeaninf-Synmean(i-1))*exp(-dt/tau_syn);     % and update synapse
end

% Now carry out the simulation with real spikes

spikes=zeros(size(t));                  % initialize spike array

for i = 1:length(t)                     % for every time point
    if rand(1) < rate(i)*dt             % probability of spike it rate(i)*dt
        spikes(i) = 1;                  % put spikes at randomly chose time points
    end
end

D = ones(size(t));                      % Depression varies with time
F = ones(size(t));                      % Facilitation varies with time
pr=p0*F;                                % Vesicle release probability increases with F
p=pr.*D;                                % Total release proportion, p, depends on D 

Ssyn = zeros(size(t));                  % initialize synaptic output array
nwidth = round(10*tau_syn1/dt); % just add current from a spike for this number of time points

spiketimes = dt*find(spikes);           % exact time of each spike
Nspikes = length(spiketimes);           % total number of spikes
for n = 1:Nspikes                       % for each spike
    ispike = round(spiketimes(n)/dt);   % find time index when spike occurs
    pr(ispike) = p0*F(ispike);          % release probability of each vesicle
    p(ispike) = pr(ispike)*D(ispike);   % total proportion of releases out of maximum possible
    
    D(ispike+1) = D(ispike)*(1-pr(ispike));     % decrease D by a factor given by proportion of releases
    pr(ispike+1) = pr(ispike)+f_f*(1-pr(ispike));    % increase probability of release by facilitation
    F(ispike+1) = pr(ispike+1)/p0;                   % update facilitation factor, F    
    p(ispike+1) = pr(ispike+1)*D(ispike+1);     % update release proportion for next time point

    if n < Nspikes                          
        nspike = round(spiketimes(n+1)/dt);     % update variables until the next spike 
    else                                        % or
        nspike = length(t);                     % until the end of the simulation
    end
    
    for i = ispike+1:nspike-1
        D(i+1) = 1 + (D(i)-1)*exp(-dt/tau_d);   % Depression decays exponentially up to 1
        F(i+1) = 1 + (F(i)-1)*exp(-dt/tau_f);   % Facilitation decays exponentially down to 1
        pr(i+1) = p0*F(i+1);                    % Update vesicle release probability
        p(i+1) = p0*D(i+1)*F(i+1);              % Update proportion of total possible synaptic conductance
    end 
    
    iend = min(ispike+nwidth,length(t))         % number of time points to update synaptic conductance
    for i=ispike+1:iend                     
        tdiff = dt*(i-ispike);                  % time since the spike
        Ssyn(i) = Ssyn(i) + ...                 % increase the synaptic conductance with saturation
            p(ispike)*(1-Ssyn(ispike))*(exp(-tdiff/tau_syn1)-exp(-tdiff/tau_syn2));        
    end
end

figure(1)
subplot(4,1,1)
plot(t,spikes)
subplot(4,1,2)
plot(t,D)
hold on
plot(t,Dmean,'r')
subplot(4,1,3)
plot(t,F)
subplot(4,1,4)
plot(t,Ssyn)
hold on
plot(t,Synmean,'r')