function [f,Conductances] =  EulerMaruyama(x,dt,ModelSwitch,NoiseSwitches,Params,NumChannelTypes,ActivationVarsPerChannelType,NumActivationVars,StatesPerChannelType,InputCurrent,NumChannelsEachType)

%get activation and inactivation rates (Alphas and Betas) for current membrane potential
%also find the A matrices (returned in a cell array, where each cell is a different channel type type)
%also find the state probabilities (returns in a cell array, p, where each cell is a different channel type)
[A,p,Alphas,Betas] = UpdateDriftAndOccupancies(x,NumActivationVars,StatesPerChannelType,ModelSwitch);

%extract the conductances from the last element of the p vectors
Conductances = zeros(1,NumChannelTypes);
for i  = 1:NumChannelTypes
    Conductances(1,i) = p{i}(end);
end

%update V and activation variables, using their previous values and fluctuations at previous time step
%for these variables, we just use the Euler method
f = x + dt*UpdateEquations(x,p,Params,NumChannelTypes,ActivationVarsPerChannelType,NumActivationVars,StatesPerChannelType,InputCurrent,Alphas,Betas);

%update the associated SDEs using the Euler-Maruyama method
j = NumActivationVars + 1;
for i  = 1:NumChannelTypes
    k = StatesPerChannelType(i);
    if NoiseSwitches(i) == 1
        %DiffusionFunc finds the necessary matrix square roots of the diffusion matrices formed from the A matrices
        f(j+1:j+k) = f(j+1:j+k) + dt * A{i}*x(j+1:j+k) + sqrt(dt) * UpdateDiffusionMatrixSquareRoot(A{i},p{i},NumChannelsEachType(i)) * randn(k,1);
    end
    j = j + k;
end