function [signals] = wc_coupled_stochastic(G,D,time,dt,c5,c6,stim_P,stim_Q)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation of coupled and stimulated Wilson-Cowan model with noise.
% Stochastic Euler-Maruyama method
% Created by Kanika Bansal and Sarah Muldoon.
% Update: Sept, 11 2018. 
% Ref: Bansal et al. Cognitive chimera states in human brain networks, Science Advance, 2019
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Description of variables
%Input:
%           'G'           --- NxN matrix of structural connectivity - note
%                             that all diagonal entries will be set to 0
%           'D'           --- NxN matrix of time delays between oscillators%                             
%           'time'        --- Total duration of simulated dynamics 
%           'dt'          --- Simulaiton step size
%           'c5'          --- The global coupling strength for excitatory population
%           'c6'          --- The global coupling strength for inhibitory populatioin
%           'stim_P'      --- nx4 matrix to define the external stimulation function to the
%                             excitatory populations where n is the number of regions (oscillators) to be stimulated. 
%                             Second dimension represents [oscillator, stimulation value, start time, end time]
%           
%           'stim_Q'      --- nx4 matrix to define the external stimulation function to the
%                             inhibitory populations where n is the number of regions (oscillators) to be stimulated. 
%                             Second dimension represents [oscillator, stimulation value, start time, end time]

%Output: 
%           signals.e     --- time x N matrix of the excitatory data, where N
%                             is the number of oscillators
%           signals.i     --- time x N matrix of the inhibitory data, where N
%                             is the number of oscillators
%           signals.t     --- vector of time 
%  

%% Equation and simulation parameters: 
N = size(G,2);                      % --- Number of brain regions (oscillators)
noise_var = 0.00005;                % --- varaiance of noise (gaussian)to be added to the oscillators

% Computational integrator used in this code is only valid for
% small noise_var.


initEstate = 0.0*randn(1,N)+0.1;    % --- size 1xN, gives the inital states of the 
                                          %excitatory population for each oscillator
initIstate = 0.0*randn(1,N)+0.1;    % --- size 1xN, gives the initial states of the
                                          %inhibitory population for each oscillator

% Model parameters obtained from Wilson HR, Cowan JD. Biophys J. 1972;12(1):1-24
tau = 8;            % --- the excitatory time constant in ms
r_e = 1;            % --- the excitatory refractory period
r_i = 1;            % --- the inhibitory refractory period
c1 =  16;           % --- parameter that defines the self-coupling of the
                    %     excitatory population
c2 = 12;            % --- parameter that defines the coupling of the
                          %inhibitory to excitatory population for a single 
                          %column
c3 = 15;            % --- parameter that defines the coupling of the
                          %excitatory to inhibitory population for a single 
                          %column
c4 = 3;             % --- parameter that defines the self-coupling of the
                          %inhibitory population
a_e = 1.3;          % --- parameter for excitatory sigmoidal
a_i = 2;            % --- parameter for inhibitory sigmoidal
theta_e = 4;        % --- parameter for excitatory sigmoidal
theta_i = 3.7;      % --- parameter for inhibitory sigmoidal


%% Setting up delay matrix, converting time delay matrix into index
delays_m_steps=round(D./dt); 


%% Simulation parameters

%setting the self_connectivity to 0 
for i=1:N
    G(i,i)=0;
end

% total simulation steps
T = (time/dt)+1;
signals.t = transpose(0:dt:time);

%parameters for sigmoid function
Smax_e=1-1/(1+exp(a_e*theta_e));
Smax_i=1-1/(1+exp(a_i*theta_i));
sig_params_e=[a_e theta_e];
sig_params_i=[a_i theta_i];


% initializing output (setting all the initial values same) 

estate = zeros(T, N);
istate= zeros(T, N);

estate(1,:) = initEstate;
istate(1,:) = initIstate;



%% Setting up stimulation matrix

%Excitatory (P)
P = zeros(1,N);          
P_matrix=repmat(P,T,1);
num_stim=size(stim_P);
num_stim=num_stim(1);
if num_stim > 0
    for i=1:num_stim
        oscillator_i=stim_P(i,1);
        amp_stim_i=stim_P(i,2);
        start_stim_i=stim_P(i,3)/dt;
        end_stim_i=stim_P(i,4)/dt;
        P_matrix(start_stim_i:end_stim_i,oscillator_i)=P_matrix(start_stim_i:end_stim_i,oscillator_i)+amp_stim_i;
    end
end

%Inhibitory (Q)
Q = zeros(1,N);
Q_matrix=repmat(Q,T,1);
num_stimQ=size(stim_Q);
num_stimQ=num_stimQ(1);
if num_stimQ > 0
    for i=1:num_stimQ
        oscillator_i=stim_Q(i,1);
        amp_stim_i=stim_Q(i,2);
        start_stim_i=stim_Q(i,3)/dt;
        end_stim_i=stim_Q(i,4)/dt;
        Q_matrix(start_stim_i:end_stim_i,oscillator_i)=Q_matrix(start_stim_i:end_stim_i,oscillator_i)+amp_stim_i;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Main Loop: Simulating coupled and stimulated WC oscillators 
%  Alert! This computational scheme is valid for small noise values.  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for t= 2:T
    %creating noise that is fed to the excitatory and inhibitory populations
    noise_e=noise_var * randn(1,N);
    noise_i=noise_var * randn(1,N);
    
    %calculating the delay time
    delayed_time_index=t-delays_m_steps;
    negative_index=find(delayed_time_index < 2);
    if ~isempty(negative_index)
        delayed_time_index(negative_index)=2;
    end
    
    %calcuating the estate and istate of the delay times
    edelayed=zeros(N);
    idelayed=zeros(N);
    for i=1:N
        edelayed_i=estate(delayed_time_index(:,i)-1,i);
        edelayed(i,:)=edelayed_i;
        idelayed_i=istate(delayed_time_index(:,i)-1,i);
        idelayed(i,:)=idelayed_i;
        
    end
   
    
    %calculating the coupling term input
    c_in_e = sum(G.*edelayed,1);
    c_in_i = sum(G.*idelayed,1);
    
    
       
    %% Calculating the states using Euler-Maruyama method 

    estate(t,:) = estate(t-1,:) + dt .* (1/tau) .* (- estate(t-1, :) + ...
                (Smax_e-r_e.*estate(t-1,:)).*sigmoidal_wc72(c1.*estate(t-1,:)-c2.*istate(t-1,:)+c5.*c_in_e+P_matrix(t,:), sig_params_e))...
                + ((1/tau).*dt^1/2.*noise_e);

    istate(t,:) = istate(t-1,:) + dt .* (1/tau) .* (- istate(t-1, :)+ ...
                (Smax_i-r_i.*istate(t-1,:)).*sigmoidal_wc72(c3.*estate(t-1,:)-c4.*istate(t-1,:)+c6.*c_in_i+Q_matrix(t,:), sig_params_i))...
                + ((1/tau).*dt^1/2.*noise_i);
        
end
signals.e = estate; 
signals.i = istate;


function [sig_val] = sigmoidal_wc72( x, sig_params )
%this is to calculate the sigmoidal funciton used in the Wilson-Cowan
%model described in Wilson and Cowan's 1972 paper

%Input:
%   'x'             --- 1xN vector where N is the
%                       number of oscillators
%   'sig_params'    --- 1x2 vector listing the 'a','theta' parameter
%                       values as in the paper

%Output: 
%   'sig_val'       --- the value of the function

a=sig_params(1);
theta=sig_params(2);

sig_val = 1./(1 + exp(- a.*(x - theta))) - 1/(1 + exp(a*theta));
