//Example computer code developed in the Scilab environment as used in:
//"Homomeric Q/R edited AMPA receptors conduct when desensitized"
//Coombs, Soto, McGee, Gold, Farrant, Cull-Candy
//Nature Communications manuscript NCOMMS-19-09071
//
//The following is code for calculating expected current waveforms and NSFA using "Scheme 1"
//The code allows reproduction of pooled data and simulated curves displayed in Figure 5 of the manuscript

//Scheme 1 [state (matrix position)]:
//         O1(6)    O2(7)     O3(8)     O4(9)
//  R(1)   R1(2)    R2(3)     R3(4)     R4(5)
//  D0(10) D1(11)   D2(12)    D3(13)    D4(14)
//                  2D2(15)   2D3(16)  2D4(17)

//Units of rates are ms^-1 except for k1on which is M^-1 ms^-1 (at 10 mM glutamate)
k1on=13      
k1off=0.35   
bet=1        
alph=3.1     
gam=0.088    
delt=0.11    
gam2=0.036    
delt2=0.039    
delt0=0.00048
gam0=0.008   
delt3=0.11   
gam3=0.69    
//single channel currents (pA at -60 mV). Divide by 0.06 to yield conductance in pS.
CondO1=0.059   
CondD1=0.01  

//Constrained by microscopic reversibility
k2off = (4*delt0*gam*k1off*3*k1on)/(4*k1on*delt3*gam0)

desnoisex=zeros(10,1)
desnoisey=zeros(10,1)
deactnoisex=zeros(10,1)
deactnoisey=zeros(10,1)
actnoisex=zeros(10,1)
actnoisex=zeros(10,1)
Importdestime=zeros(27,1)
Desenscurrent=zeros(27,1)
Importdeacttime=zeros(25,1)
Deactcurrent=zeros(25,1)
desnoisepoints=zeros(10,1)
deactnoisepoints=zeros(10,1)
Importacttime=zeros(10,1)
Actcurrent=zeros(10,1)

//Global averaged parsed desensitizing current.
Importdestime = [-5,0,0.4,0.7,0.9,1,1.7,2.6,3.7,4.9,6.7,9.1,13.7,21.7,44.7,52.7,62.7,72.7,82.7,100,100.6,101,101.8,102.9,103.5,104.5,106,107.5,111,114,118,124,134,156,176]
Desenscurrent = [0,0,-0.5,-1,-0.99,-0.978595,-0.91561,-0.863629,-0.807657,-0.754217,-0.697109,-0.643016,-0.584499,-0.531182,-0.473596,-0.466,-0.4605,-0.457,-0.4555,-0.4538,-0.43,-0.4,-0.35,-0.28,-0.25,-0.22,-0.19,-0.16,-0.13,-0.1,-0.08,-0.06,-0.04,-0.02,-0.012]

//Global averaged parsed deactivating current.
Importdeacttime = [-3,0,0.05,0.1,0.14,0.16,0.18,0.21,0.25,0.29,0.37,0.5,0.65,0.75,1.1,1.5,2,2.4,2.9,3.6,4.7,7.6,9.6,11.8,16.1,21.1,30.1,40.1,47.4]
Deactcurrent = [0,0,-0.124669,-0.222272,-0.328289,-0.426022,-0.528783,-0.630926,-0.731496,-0.830402,-0.932031,-1,-0.96,-0.935026,-0.831095,-0.733975,-0.632953,-0.535268,-0.428732,-0.326188,-0.221909,-0.118769,-0.09,-0.07,-0.05,-0.03,-0.0160628,-0.01,-0.005]

//Activation current (taken from deactivating data)
Importacttime = [-0.1,0,0.05,0.1,0.14,0.16,0.18,0.21,0.25,0.29,0.37,0.5,0.65] 
Actcurrent = [0,0,-0.124669,-0.222272,-0.328289,-0.426022,-0.528783,-0.630926,-0.731496,-0.830402,-0.932031,-1,-0.96]

//Current(x)-variance(y) relationships (baseline subtracted) for the above currents
desnoisex = [0.978595,0.91561,0.863629,0.807657,0.754217,0.697109,0.643016,0.584499,0.531182,0.473597]
desnoisey = [0.145556,0.139898,0.124039,0.111702,0.096716,0.086429,0.070510,0.054737,0.042748,0.033389]
deactnoisex = [0.93203,0.83040,0.73150,0.63093,0.52878,0.42602,0.32829,0.22227,0.12467,0.02029]
deactnoisey = [0.12302,0.12272,0.10146,0.08527,0.05368,0.03536,0.02849,0.01310,0.00784,0.00479]
actnoisex = [0.93503,0.83109,0.73397,0.63295,0.53527,0.42873,0.32619,0.22191,0.11877,0.01606]
actnoisey = [0.13152,0.11515,0.14047,0.10543,0.09344,0.08189,0.06591,0.04828,0.02832,0.00450]

subplot(231)
plot2d(Importdestime,Desenscurrent)
subplot(234)
plot2d(desnoisex,desnoisey)
subplot(232)
plot2d(Importdeacttime,Deactcurrent)
subplot(235)
plot2d(deactnoisex,deactnoisey)
subplot(233)
plot2d(Importacttime,Actcurrent)
subplot(236)
plot2d(actnoisex,actnoisey)

//This first section sets up a Q-matrix to define occupancies of the four states
Q = zeros(17,17)    //Creates a blank matrix for transitions

Q(1,2) = 0     //Q is in the absence of agonist hence this is zero
Q(1,10) = 4*delt0
Q(2,1) = k1off
Q(2,3) = 0
Q(2,6) = bet
Q(2,11) = delt3  
Q(3,2) = 2*k1off  
Q(3,4) = 0  
Q(3,7) = 2*bet  
Q(3,12) = 2*delt 
Q(4,3) = 3*k1off
Q(4,5) = 0
Q(4,8) = 3*bet
Q(4,13) = 3*delt
Q(5,4) = 4*k1off
Q(5,9) = 4*bet
Q(5,14) = 4*delt
Q(6,2) = alph
Q(7,3) = alph
Q(7,6) = k1off
Q(8,4) = alph
Q(8,7) = 2*k1off
Q(9,5) = alph
Q(9,8) = 3*k1off
Q(10,1) = gam0
Q(11,2) = gam
Q(11,10) = k2off
Q(11,12) = 0
Q(12,3) =gam
Q(12,11) =k1off
Q(12,13) =0
Q(12,15) =delt2
Q(13,4) =gam
Q(13,12) =2*k1off
Q(13,14) = 0
Q(13,16) =2*delt2
Q(14,5) =gam
Q(14,13) =3*k1off
Q(14,17) =3*delt2
Q(15,12) =gam2
Q(15,16) = 0
Q(16,13) =gam2
Q(16,15) =k1off
Q(16,17) = 0
Q(17,14) =gam2
Q(17,16) =2*k1off

//the sum of each row needs to be zero for the Q-matrix so the diagonals are:
Q(1,1) = -Q(1,2)-Q(1,10)
Q(2,2) = -Q(2,1)-Q(2,3)-Q(2,6)-Q(2,11)
Q(3,3) = -Q(3,2)-Q(3,4)-Q(3,7)-Q(3,12)
Q(4,4) = -Q(4,3)-Q(4,5)-Q(4,8)-Q(4,13)
Q(5,5) = -Q(5,4)-Q(5,9)-Q(5,14)
Q(6,6) = -Q(6,2)
Q(7,7) = -Q(7,3)-Q(7,6)
Q(8,8) = -Q(8,4)-Q(8,7)
Q(9,9) = -Q(9,5)-Q(9,8)
Q(10,10) = -Q(10,1)-Q(10,11)
Q(11,11) = -Q(11,2)-Q(11,10)-Q(11,12)
Q(12,12) = -Q(12,3)-Q(12,11)-Q(12,13)-Q(12,15)
Q(13,13) = -Q(13,4)-Q(13,12)-Q(13,14)-Q(13,16)
Q(14,14) = -Q(14,5)-Q(14,13)-Q(14,17)
Q(15,15) = -Q(15,12)-Q(15,16)
Q(16,16) = -Q(16,13)-Q(16,15)-Q(16,17)
Q(17,17) = -Q(17,14)-Q(17,16)

Q2 = zeros(17,17)  //Q2 is the Q-matrix in the presence of agonist 
Q2 = Q           //Mostly identical to Q
Q2(1,2) = 4*k1on                     
Q2(1,1) = Q(1,1)-Q2(1,2)
Q2(2,3) = 3*k1on
Q2(2,2) = Q(2,2)-Q2(2,3)
Q2(3,4) = 2*k1on
Q2(3,3) = Q(3,3)-Q2(3,4)
Q2(4,5) = k1on
Q2(4,4) = Q(4,4)-Q2(4,5)
Q2(6,7) = 3*k1on
Q2(6,6) = Q(6,6)-Q2(6,7)
Q2(7,8) = 2*k1on
Q2(7,7) = Q(7,7)-Q2(7,8)
Q2(8,9) = k1on
Q2(8,8) = Q(8,8)-Q2(8,9)
Q2(10,11) = 3*k1on 
Q2(10,10) = Q(10,10)-Q2(10,11)
Q2(11,12) = 3*k1on 
Q2(11,11) = Q(11,11)-Q2(11,12)
Q2(12,13) = 2*k1on 
Q2(12,12) = Q(12,12)-Q2(12,13)
Q2(13,14) = k1on 
Q2(13,13) = Q(13,13)-Q2(13,14)
Q2(15,16) = 2*k1on 
Q2(15,15) = Q(15,15)-Q2(15,16)
Q2(16,17) = k1on 
Q2(16,16) = Q(16,16)-Q2(16,17)

p = zeros(1,17)

I = eye(17,18) 
S = zeros(17,18)
S = Q*I
S(:,$) = 1
St = S'
u = zeros(1,17)
u = [ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
p = zeros(1,17)
p = u/(S*St)  //Initial occupancies

//Calculating occupancies for the desensitization  pulse
y = zeros(35,18)                        //Creates a table, y, to store occupanicies for each state at time t

for t = [3:1:20]                        //A loop for each timepoint of glutamate application
    x = [p*expm(Q2*Importdestime(t))],  //Calculates matrix with occupancies (x) at timepoint from Importdestime
    for z = [2:1:18]                    //This loop puts the occupancies of each state (z) at time t, in table y
        y(t,z) = x(z-1)
        y(t,1) = Importdestime(t)       //lists time point in column 1 of the matrix
    end    
end  

a = p*expm(Q2*Importdestime(t))  //occupancies at end of previous step into 'a'

//The following works out the occupancies following glutamate removal

for t = [21:1:35]
    x =[a*expm(Q*(Importdestime(t)-(Importdestime(20))))],  //The time since glutamate removal. 
    for z = [2:1:18]                                        //as above
        y(t,z) = x(z-1) 
        y(t,1) = Importdestime(t) 
    end 
end                              //takes into account 20ms before application

///////////////////////////////////////////////////////////////////
//Expected current from above occupancies
O = zeros(35,1)
D = zeros(35,1)
O2 = zeros(35,1)
D2 = zeros(35,1)
O3 = zeros(35,1)
D3 = zeros(35,1)
O4 = zeros(35,1)
D4 = zeros(35,1)

for t = [1:1:35]
O(t) = y(t,7)           //State O1
D(t) = y(t,12)    //State D1
O2(t) = y(t,8)         //State O2
D2(t) = y(t,13)+y(t,16) //States D2, 2D2
O3(t) = y(t,9)         //State O3
D3(t) = y(t,14)+y(t,17) //States D3, 2D3
O4(t) = y(t,10)         //State O4
D4(t) = y(t,15)+y(t,18) //States D4, 2D4
end

plty = zeros(35,1) //For the current at each time

//Calculating the number of channels given the peak current, occupanices and unitary currents
Ndes=-Desenscurrent(5)/(CondO1*O(5)+CondD1*D(5)+2*CondO1*O2(5)+2*CondD1*D2(5)+3*CondO1*O3(5)+3*CondD1*D3(5)+4*CondO1*O4(5)+4*CondD1*D4(5))

for t = [1:1:35]  
    plty(t,1) = -Ndes*(CondO1*O(t)+CondD1*D(t)+2*CondO1*O2(t)+2*CondD1*D2(t)+3*CondO1*O3(t)+3*CondD1*D3(t)+4*CondO1*O4(t)+4*CondD1*D4(t))
end  

//Initialize noise matricies
Noisex = zeros (1,5) 
Noisey = zeros (1,5) 

for t = [6:1:15]
//The noise is expected to be N(po)i^2(1-po) for each occupancy, and the total noise is the sum for each state
    Noisey(1,t-5)= Ndes*(O(t)*CondO1*CondO1*(1-O(t))+D(t)*CondD1*CondD1*(1-D(t))+O2(t)*2*CondO1*2*CondO1*(1-O2(t))+D2(t)*2*CondD1*2*CondD1*(1-D2(t))+O3(t)*9*CondO1*CondO1*(1-O3(t))+D3(t)*9*CondD1*CondD1*(1-D3(t))+O4(t)*16*CondO1*CondO1*(1-O4(t))+D4(t)*16*CondD1*CondD1*(1-D4(t)))
    Noisex(1,t-5)=-plty(t,1)
end  

///////////////////////////////////////////////////////////////////
//Deactivation occupancies

y2 = zeros(29,18)  //Creates a table, y2, to store deactivation/activation occupanicies at time t

//Calculates occupancies during activation
for t = [3:1:15] 
  x = [p*expm(Q2*Importdeacttime(t))], //Calculating occupancies following the jump into glutamate
  for z = [2:1:18]      //As for desensitization above
    y2(t,z) = x(z-1)
    y2(t,1) = Importdeacttime(t) 
  end    //lists time point in column 1 of the mat
end

a = [p*expm(Q2*(Importdeacttime(t)))] //occupancies at end of previous step into 'a'

//Calculates occupancies during deactivation
for t = [16:1:29]
    x =[a*expm(Q*(Importdeacttime(t)-(Importdeacttime(15))))],  
    for z = [2:1:18]      
        y2(t,z) = x(z-1) 
        y2(t,1) = Importdeacttime(t) 
    end 
end  

///////////////////////////////////////////////////////////////////
//Expected current from above occupancies
Oa = zeros(29,1)
Da = zeros(29,1)
O2a = zeros (29,1)
D2a = zeros (29,1)

for t = [1:1:29]
Oa(t) = y2(t,7)           //State O
Da(t) = y2(t,12)    //States D1
O2a(t) = y2(t,8)         //State O2
D2a(t) = y2(t,13)+y2(t,16) //States D2, 2D2
O3a(t) = y2(t,9)         //State O2
D3a(t) = y2(t,14)+y2(t,17) //States D3, 3D2
O4a(t) = y2(t,10)         //State O2
D4a(t) = y2(t,15)+y2(t,18) //States D4, 4D2
end

plty2 = zeros(29,1) //Deact current
plty3 = zeros(13,1) //Act current

//Calculating the number of channels given the peak current, occupanices and unitary currents
Ndeact=((-Deactcurrent(12)/(CondO1*Oa(12)+CondD1*Da(12)+2*CondO1*O2a(12)+2*CondD1*D2a(12)+3*CondO1*O3a(12)+3*CondD1*D3a(12)+4*CondO1*O4a(12)+4*CondD1*D4a(12)))+(-Deactcurrent(13)/(CondO1*Oa(13)+CondD1*Da(13)+2*CondO1*O2a(13)+2*CondD1*D2a(13)+3*CondO1*O3a(13)+3*CondD1*D3a(13)+4*CondO1*O4a(13)+4*CondD1*D4a(13))))/2

//Calculating the current
for t = [1:1:29]
    plty2(t,1) = -Ndeact*(CondO1*Oa(t)+CondD1*Da(t)+2*CondO1*O2a(t)+2*CondD1*D2a(t)+3*CondO1*O3a(t)+3*CondD1*D3a(t)+4*CondO1*O4a(t)+4*CondD1*D4a(t))
end  
for t=[1:1:13]
    plty3(t,1) = plty2(t,1)
end

Noise2y = zeros (1,10) 
Noise3y = zeros (1,10) 
Noise2x = zeros (1,10) 
Noise3x = zeros (1,10) 

//The noise is expected to be N(po)i^2(1-po) for each occupancy, and the total noise is the sum for each state
//activation noise is calculated from current timepoints 3-12
for t = [3:1:12] 
Noise3y(1,t-1)= Ndeact*(Oa(t)*CondO1*CondO1*(1-Oa(t))+Da(t)*CondD1*CondD1*(1-Da(t))+O2a(t)*2*CondO1*2*CondO1*(1-O2a(t))+D2a(t)*2*CondD1*2*CondD1*(1-D2a(t))+O3a(t)*9*CondO1*CondO1*(1-O3a(t))+D3a(t)*9*CondD1*CondD1*(1-D3a(t))+O4a(t)*16*CondO1*CondO1*(1-O4a(t))+D4a(t)*16*CondD1*CondD1*(1-D4a(t)))
Noise3x(1,t-1)=-plty2(t,1)
end

//activation noise is calculated from current timepoints 14-22 and 27
for t = [14:1:22]
Noise2y(1,t-13)= Ndeact*(Oa(t)*CondO1*CondO1*(1-Oa(t))+Da(t)*CondD1*CondD1*(1-Da(t))+O2a(t)*2*CondO1*2*CondO1*(1-O2a(t))+D2a(t)*2*CondD1*2*CondD1*(1-D2a(t))+O3a(t)*9*CondO1*CondO1*(1-O3a(t))+D3a(t)*9*CondD1*CondD1*(1-D3a(t))+O4a(t)*16*CondO1*CondO1*(1-O4a(t))+D4a(t)*16*CondD1*CondD1*(1-D4a(t)))
Noise2x(1,t-13)=-plty2(t,1)
end

t=27
Noise2y(1,10)= Ndeact*(Oa(t)*CondO1*CondO1*(1-Oa(t))+Da(t)*CondD1*CondD1*(1-Da(t))+O2a(t)*2*CondO1*2*CondO1*(1-O2a(t))+D2a(t)*2*CondD1*2*CondD1*(1-D2a(t))+O3a(t)*9*CondO1*CondO1*(1-O3a(t))+D3a(t)*9*CondD1*CondD1*(1-D3a(t))+O4a(t)*16*CondO1*CondO1*(1-O4a(t))+D4a(t)*16*CondD1*CondD1*(1-D4a(t)))
Noise2x(1,10)=-plty2(t,1)

subplot(231) //plots the top left of 4x2
plot2d (Importdestime,plty,5,leg="desensitization") //plotting the desensitizing response

subplot(232) //plots the top left of 4x2
plot2d (Importdeacttime,plty2,5,leg="deactivation") //plotting the desensitizing response

subplot(233) //plots the top left of 4x2
plot2d (Importacttime,plty3,5,leg="activation") //plotting the desensitizing response


subplot(234) //plots noise under desensitization current plot
plot2d(Noisex,Noisey,5,leg="desensitization noise");

a=gca() //Set axes from 0 to max to include origin.
a.data_bounds=[0,0;max(Noisex),max(Noisey)] 

subplot(235) //plots noise under desensitization current plot
plot2d(Noise2x,Noise2y,5,leg="deactivation noise");

subplot(236) //plots noise under desensitization current plot
plot2d(Noise3x,Noise3y,5,leg="activation noise");


