function [curr_voltage,AP]=LIF_one_step_multiple_neurons(curr_voltage,Input,params)
%LIF_ONE_STEP_MULTIPLE_NEURONS Runs one step of simulation - LIF neurons
%% Model Parameters
if nargin<3
    params.tau=15;                 % time constant (ms)
    params.input_gain=100;         % gain to scale input (normalized xcorr is [-1,1])
    params.noise_sd=0;             % number of SDs for gaussian input noise
    params.AP_threshold=-30;       % AP threshold (mV)
    params.num_steps_per_ms=fs_ms; % num of steps per ms
    params.lateral_inhibition=5;   % gain to scale lateral inhibition
    params.rest_potential=-70;     % resting state potential (mV)
end

%% Set important variables
I_ext=params.input_gain .* ... scale ->
    (Input+randn(size(Input))*params.noise_sd); %input + noise
tau=params.tau*params.num_steps_per_ms; %update tau constant for step size

AP=zeros(size(Input));
%% One Step Integration with Euler method

curr_firing=curr_voltage>params.AP_threshold; %which neurons currently fire?
curr_voltage(curr_firing)=params.rest_potential; %reset voltage
AP(curr_firing)=1; %mark AP at for those that fire

%not firing - update voltage:
curr_voltage(~curr_firing)=curr_voltage(~curr_firing)+ ... %the current voltage
    (params.rest_potential-curr_voltage(~curr_firing))./tau+... %minus leak
    I_ext(~curr_firing)./tau -... %plus (noisy?) input
    sum(AP)*params.lateral_inhibition/tau; %minus lateral inhibition

end

