function maxima_val=single_syn_stsd(AP_times)
%the destiny of this m-file is making predictions about [Ca]bouton, exocytosis of transmitter, [glutamate]cleft, and AMPAR conductance(s), among other things
%close all
plot_it=1

tic
%global waveform_holder     %used for collecting data from multiple trials
%waveform_holder=[]

%AP_times=[200 400]                      %deterministic train of APs
% AP_times=ceil(cumsum(exprnd(100,1,3)))     %stochastic train of APs
trial_duration=1000;    %a timestep is roughly a millisecond

calcium_conc=.005;      %in presynaptic cytoplasm (starting with above zero makes spontaneous release possible)
calcium_persistence=0.97;%the factor by which current value is multiplied each time step (larger values, slower decay)
RRP_size_max=4;        %in vesicles
RRP_size=RRP_size_max;
RP_to_RRP_rate=.003;  %in vesicles per time step
RP_size=100;            %in vesicles
PE_size=0;              %in vesicles (the post-exocytotic pool)
PE_to_RP_rate=.005;       %in vesicles per time step
PE_to_RRP_rate=0;       %in vesicles per time step %RP to RRP and PE to RRP sum to the RRP refilling rate
glu_per_vesicle=5;      %in effective pA(?)
glu_cleft=.005;         %in effective pA(?)
glutamate_cleft_persistence=.8;
exocytosis_unlikelihood=10;  %this is multiplied by a random number -- if the result is greater than [Ca]bouton, release is impossible

number_of_ampars=3000;
ampar_bound=zeros(1,number_of_ampars);  %at the beginning, none is bound
ampar_open=ampar_bound;                 %therefore none is open
ampar_conductance=ampar_bound;          %therefore there is no conductance
initial_ampar_conductance=.1;   %in pA/time step
ampar_open_threshold=.7;    %used for determining whether a bound, closed ampar moves to the bound, open state
ampar_unbind_threshold=.3;  %used for determining whether a bound, closed ampar moves to the unbound, closed state
ampar_close_threshold=.1;   %used for determining whether a bound,open ampar moves to the bound, closed state
ampar_binding_unlikelihood=20; %larger values make it less likely that an unbound receptor will bind glutamate at a given [glu]cleft
ampar_conductance_persistence=.9;
synaptic_conductance=0;

calcium_hx=[];
RRP_size_hx=[];
RP_to_RRP_rate_hx=[];
RP_size_hx=[];
PE_size_hx=[];
PE_to_RP_rate_hx=[];
PE_to_RRP_rate_hx=[];
glu_cleft_hx=[];
synaptic_conductance_hx=[];

exocytosis_permitted=1;
for i=1:trial_duration
    
    %APs drive [calcium] in the bouton (aka the release probability) to 1
    if find(i==AP_times)
        calcium_conc=calcium_conc+1     %summation across APs can lead to facilitation
    else calcium_conc=calcium_conc*calcium_persistence;
    end
    
    if calcium_conc<.1, exocytosis_permitted=1; end %this serves to make release, if prevented immediately after release, possible again, once the [Ca] falls low enough (i.e. if ample time has passed)
    %stochastic exocytosis based on the calcium concentration, as long as some other criteria are met
    if (calcium_conc-(exocytosis_unlikelihood*rand))>0 & exocytosis_permitted & RRP_size>=1 %a release event    -- should release affect the probability of release? (uniquantal rule) --> c.f. if calcium_conc<n, exocytosis_permitted=1; end (above)
        exocytosis_permitted=1;     %set to zero to prevent additional release until the calcium drops *low* enough (surrogate for how much time has passed since the release)
        RRP_size=RRP_size-1;
        glu_cleft=glu_cleft+glu_per_vesicle;
        PE_size=PE_size+1;           %(increments each time there is a release, until retrieved) -- should there be kiss-and-run (exocytosis w/o vesicle fusion), this best not increment
    end
    
    %the effect of the cleft glutamate on ampa receptor binding and opening
    for j=1:number_of_ampars
        if ampar_bound(j)==0
            if ampar_open(j),error('ydfu'),end     %a check to make sure unbound receptors are not opening
            if (glu_cleft/(ampar_binding_unlikelihood*glu_per_vesicle))>rand   %the glu_per_vesicle is used as an imperfect normalizer, a proxy for max possible [glu]cleft
                ampar_bound(j)=1;   %binding glutamate
            end
		elseif ampar_bound(j)==1&ampar_open(j)==0
            randi=rand;
            if randi>ampar_open_threshold         %moves to the open state
                ampar_open(j)=1;
            elseif  randi<ampar_unbind_threshold    %moves back to unbound state
                ampar_bound(j)=0;
            end
		elseif ampar_bound(j) & ampar_open(j)
            if rand<ampar_close_threshold
                ampar_open(j)=0; ampar_conductance(j)=0; %moves back to the closed state
            elseif ampar_conductance(j)
                ampar_conductance(j)=ampar_conductance(j)*ampar_conductance_persistence;  %if the receptor is already open, it desensitizes one unit
            else
                ampar_conductance(j)=initial_ampar_conductance;   %newly opened receptors become open (conducting) on the next time step
            end
        end
    end
    synaptic_conductance=sum(ampar_conductance);
    
    %should activity change the rate of refilling?
    
    %changes performed each timestep
    if (RRP_size + RP_to_RRP_rate +PE_to_RRP_rate) <= RRP_size_max, RRP_size = RRP_size + RP_to_RRP_rate +PE_to_RRP_rate;end
    RP_size=RP_size - RP_to_RRP_rate + PE_to_RP_rate;   %no maximum size here (should this pool not decrement when the RRP is full?)
    PE_size=PE_size - PE_to_RP_rate - PE_to_RRP_rate;
    glu_cleft=glu_cleft*glutamate_cleft_persistence;
    
    %record keeping
    calcium_hx(i)=calcium_conc;
    RRP_size_hx(i)=RRP_size;
	RP_to_RRP_rate_hx(i)=RP_to_RRP_rate;
	RP_size_hx(i)=RP_size;
	PE_size_hx(i)=PE_size;
	PE_to_RP_rate_hx(i)=PE_to_RP_rate;
	PE_to_RRP_rate_hx(i)=PE_to_RRP_rate;
	glu_cleft_hx(i)=glu_cleft;   
    synaptic_conductance_hx(i)=synaptic_conductance;
end

%finding maxima
maxima_ind=[];maxima_val=[];
[maxima_val(1) maxima_ind(1)]=max(synaptic_conductance_hx(AP_times(1):AP_times(1)+45));
[maxima_val(2) maxima_ind(2)]=max(synaptic_conductance_hx(AP_times(2):AP_times(2)+45));
maxima_ind(1)=maxima_ind(1)+AP_times(1);
maxima_ind(2)=maxima_ind(2)+AP_times(2);

if plot_it
	figure(1)
	set(gcf,'position',[26   162   1269   720])
	subplot(2,2,1)
	plot(calcium_hx,'-ko'), title('[Ca]_b_o_u_t_o_n', 'fontsize',30)
	subplot(2,2,2)
	plot(glu_cleft_hx,'-rs'), title('[glu]_c_l_e_f_t','fontsize',30)
	subplot(2,2,3)
	plot(synaptic_conductance_hx,'-b^'),title('G_s_y_n_a_p_s_e','fontsize',30)
	subplot(2,2,4)
	plot(RRP_size_hx,'-gh'),title('RRP size','fontsize',30)
	subplot(2,2,3), hold on, plot(maxima_ind(1), maxima_val(1), 'o','markersize',20),plot(maxima_ind(2),maxima_val(2),'o','color','r','markersize',20)
    drawnow
end
%waveform_holder=[waveform_holder;synaptic_conductance_hx]; %used for grabbing data from multiple trials
toc