% This program is designed for calculating the free energy surface (F(x)) during the
% process of fibril aggregation (hexamer by peptide GNNQQNY in current case). The main procedure is based on the free energy fast descent method
% introduced in the Supporting Informatnion.
% This code is programed by Liu Hong, Xianghong Qi and Yang Zhang in University of Michigan in 2011.

% Miyazawa-Jernigan parameters for the strength of inter-residue interaction (1996, JMB, 256, P623-644)
% Order of residues in the matrix: 1 Cys 2 Met 3 Phe 4 Ile 5 Leu 6 Val 7
% Trp 8 Tyr 9 Ala 10 Gly 11 Thr 12 Ser 13 Asn 14 Gln 15 Asp 16 Glu 17 His 18 Arg 19 Lys 20 Pro
WJP=[5.44,	4.99,	5.8,	5.5,	5.83,	4.96,	4.95,	4.16,	3.57,	3.16,	3.11,	2.86,	2.59,	2.85,	2.41,	2.27,	3.6,	2.57,	1.95,	3.07;
     4.99,	5.46,	6.56,	6.02,	6.41,	5.32,	5.55,	4.91,	3.94,	3.39,	3.51,	3.03,	2.95,	3.3,	2.57,	2.89,	3.98,	3.12,	2.48,	3.45;
     5.8,	6.56,	7.26,	6.84,	7.28,	6.29,	6.16,	5.66,	4.81,	4.13,	4.28,	4.02,	3.75,	4.1,	3.48,	3.56,	4.77,	3.98,	3.36,	4.25;
     5.5,   6.02,	6.84,	6.54,	7.04,	6.05,	5.78,	5.25,	4.58,	3.78,	4.03,	3.52,	3.24,	3.67,	3.17,	3.27,	4.14,	3.63,	3.01,	3.76;
     5.83,	6.41,	7.28,	7.04,	7.37,	6.48,	6.14,	5.67,	4.91,	4.16,	4.34,	3.92,	3.74,	4.04,	3.4,	3.59,	4.54,	4.03,	3.37,	4.2;
     4.96,	5.32,	6.29,	6.05,	6.48,	5.52,	5.18,	4.62,	4.04,	3.38,	3.46,	3.05,	2.83,	3.07,	2.48,	2.67,	3.58,	3.07,	2.49,	3.32;
     4.95,	5.55,	6.16,	5.78,	6.14,	5.18,	5.06,	4.66,	3.82,	3.42,	3.22,	2.99,   3.07,	3.11,	2.84,	2.99,	3.98,	3.41,	2.69,	3.73;
     4.16,	4.91,	5.66,	5.25,	5.67,	4.62,	4.66,	4.17,	3.36,	3.01,	3.01,	2.78,	2.76,	2.97,	2.76,	2.79,	3.52,	3.16,	2.6,	3.19;
     3.57,	3.94,	4.81,	4.58,	4.91,	4.04,	3.82,	3.36,	2.72,	2.31,	2.32,	2.01,	1.84,	1.89,	1.7,	1.51,	2.41,	1.83,	1.31,	2.03;
     3.16,	3.39,	4.13,	3.78,	4.16,	3.38,	3.42,	3.01,	2.31,	2.24,	2.08,	1.82,	1.74,	1.66,	1.59,	1.22,	2.15,	1.72,	1.15,	1.87;
     3.11,	3.51,	4.28,	4.03,	4.34,	3.46,	3.22,	3.01,	2.32,	2.08,	2.12,	1.96,	1.88,	1.9,	1.8,	1.74,	2.42,	1.9,	1.31,	1.9;
     2.86,	3.03,	4.02,	3.52,	3.92,	3.05,	2.99,	2.78,	2.01,	1.82,	1.96,	1.67,	1.58,	1.49,	1.63,	1.48,	2.11,	1.62,	1.05,	1.57;
     2.59,	2.95,	3.75,	3.24,	3.74,	2.83,	3.07,	2.76,	1.84,	1.74,	1.88,	1.58,	1.68,	1.71,	1.68,	1.51,	2.08,	1.64,	1.21,	1.53;
     2.85,	3.3,	4.1,	3.67,	4.04,	3.07,	3.11,	2.97,	1.89,	1.66,	1.9,	1.49,	1.71,	1.54,	1.46,	1.42,	1.98,	1.8,	1.29,	1.73;
     2.41,	2.57,	3.48,	3.17,	3.4,	2.48,	2.84,	2.76,	1.7,	1.59,	1.8,	1.63,	1.68,	1.46,	1.21,	1.02,	2.32,	2.29,	1.68,	1.33;
     2.27,	2.89,	3.56,	3.27,	3.59,	2.67,	2.99,	2.79,	1.51,	1.22,	1.74,	1.48,	1.51,	1.42,	1.02,	0.91,	2.15,	2.27,	1.8,	1.26;
     3.6,	3.98,	4.77,	4.14,   4.54,	3.58,	3.98,	3.52,	2.41,	2.15,	2.42,	2.11,	2.08,	1.98,	2.32,	2.15,	3.05,	2.16,	1.35,	2.25;
     2.57,	3.12,	3.98,	3.63,	4.03,	3.07,	3.41,	3.16,	1.83,	1.72,	1.9,	1.62,	1.64,	1.8,	2.29,	2.27,	2.16,	1.55,	0.59,	1.7;
     1.95,	2.48,	3.36,	3.01,	3.37,	2.49,	2.69,	2.6,	1.31,	1.15,	1.31,	1.05,	1.21,	1.29,	1.68,	1.8,	1.35,	0.59,	0.12,	0.97;
     3.07,	3.45,	4.25,	3.76,	4.2,    3.32,	3.73,	3.19,	2.03,	1.87,	1.9,	1.57,	1.53,	1.73,	1.33,	1.26,	2.25,	1.7,    0.97,	1.75];

%unit matrix for constructing hexamer:
I6=[1, 0, 0, 0, 0, 0;
    0, 1, 0, 0, 0, 0;
    0, 0, 1, 0, 0, 0;
    0, 0, 0, 1, 0, 0;
    0, 0, 0, 0, 1, 0;
    0, 0 ,0, 0, 0, 1];
% unit matrix for a peptide with 7 residues
I7=[1, 0, 0, 0, 0, 0, 0;
    0, 1, 0, 0, 0, 0, 0;
    0, 0, 1, 0, 0, 0, 0;
    0, 0, 0, 1, 0, 0, 0;
    0, 0, 0, 0, 1, 0, 0;
    0, 0 ,0, 0, 0, 1, 0;
    0, 0, 0, 0, 0, 0, 1];
%Contact matrix for parallel peptides in hexamer:
I6p=[0, 0, 1, 0, 0, 0;
    0, 0, 0, 1, 0, 0;
    1, 0, 0, 0, 1, 0;
    0, 1, 0, 0, 0, 1;
    0, 0, 1, 0, 0, 0;
    0, 0 ,0, 1, 0, 0];
%Contact matrix for antiparallel peptides in hexamer:
I6a=[0, 1, 0, 0, 0, 0;
     1, 0, 1, 0, 0, 0;
     0, 1, 0, 1, 0, 0;
     0, 0, 1, 0, 1, 0;
     0, 0, 0, 1, 0, 1;
     0, 0 ,0, 0, 1, 0];
%Rouse matrix for a nearest neighbor harmonic chain:
K=[ 1,-1, 0, 0, 0, 0, 0;
   -1, 2,-1, 0, 0, 0, 0;
    0,-1, 2,-1, 0, 0, 0;
    0, 0,-1, 2,-1, 0, 0;
    0, 0, 0,-1, 2,-1, 0;
    0, 0 ,0, 0,-1, 2,-1;
    0, 0, 0, 0, 0,-1, 1];
%Boundary effect:
Delta=[ 1,-1, 0, 0, 0, 0, 0;
       -1, 1, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 0, 0;
        0, 0, 0, 0, 0, 1,-1;  
        0, 0, 0, 0, 0,-1, 1];
%Chain stiffness g:
g=0.8;
%Inverse of the monomer correlations Xsi:
Xsi=(1-g)/(1+g)*K+g/(1-g*g)*K*K-g*g/(1-g*g)*Delta;

%Amino acide sequence for peptide GNNQQNY in hexamer
Seq=[10, 13, 13, 14, 14, 13, 8, 10, 13, 13, 14, 14, 13, 8, 10, 13, 13, 14, 14, 13, 8, 10, 13, 13, 14, 14, 13, 8, 10, 13, 13, 14, 14, 13, 8, 10, 13, 13, 14, 14, 13, 8];
%Coordinates for C_a residues obtain from PDB files.
rn=[ -9.052,   4.207,   4.651;  
     -6.522,   2.038,   2.831;  
     -3.193,   1.904,   4.589;  
      0.384,   1.888,   3.199;  
      3.270,   2.361,   5.640;  
      6.831,   2.310,   4.318;  
      9.159,   2.144,   7.299;%% Peptide 1
      9.052,   6.640,  -4.651;  
      6.522,   4.471,  -2.831;  
      3.193,   4.337,  -4.589;  
     -0.384,   4.321,  -3.199;  
     -3.270,   4.794,  -5.640;  
     -6.831,   4.743,  -4.318;  
     -9.159,   4.577,  -7.299;%% Peptide 2
     -9.052,   9.073,   4.651; 
     -6.522,   6.904,   2.831;  
     -3.193,   6.770,   4.589;  
      0.384,   6.754,   3.199; 
      3.270,   7.227,   5.640;  
      6.831,   7.176,   4.318;  
      9.159,   7.010,   7.299;%% Peptide 3
      9.052,  11.506,  -4.651;  
      6.522,   9.337,  -2.831;  
      3.193,   9.203,  -4.589;  
     -0.384,   9.187,  -3.199;  
     -3.270,   9.660,  -5.640;  
     -6.831,   9.609,  -4.318;  
     -9.159,   9.443,  -7.299;%% Peptide 4     
     -9.052,  13.939,   4.651;  
     -6.522,  11.770,   2.831;  
     -3.193,  11.636,   4.589;  
      0.384,  11.620,   3.199;  
      3.270,  12.093,   5.640;  
      6.831,  12.042,   4.318;  
      9.159,  11.876,   7.299;%% Peptide 5        
      9.052,  16.372,  -4.651;  
      6.522,  14.203,  -2.831;  
      3.193,  14.069,  -4.589;  
     -0.384,  14.053,  -3.199;  
     -3.270,  14.526,  -5.640;  
     -6.831,  14.475,  -4.318;  
     -9.159,  14.309,  -7.299%% Peptide 6
     ];

% intra-peptide interaction between residues i and j(j>=i+3) 
M1=[0, 0, 0, 1, 1, 1, 1;
    0, 0, 0, 0, 1, 1, 1;
    0, 0, 0, 0, 0, 1, 1;
    1, 0, 0, 0, 0, 0, 1;
    1, 1, 0, 0, 0, 0, 0;
    1, 1 ,1, 0, 0, 0, 0;
    1, 1, 1, 1, 0, 0, 0];
% inter-peptide interaction between residues i and j (|j|<=i+1) in parallel peptides:
M3=[1, 1, 0, 0, 0, 0, 0;
    1, 1, 1, 0, 0, 0, 0;
    0, 1, 1, 1, 0, 0, 0;
    0, 0, 1, 1, 1, 0, 0;
    0, 0, 0, 1, 1, 1, 0;
    0, 0 ,0, 0, 1, 1, 1;
    0, 0, 0, 0, 0, 1, 1];
% inter-peptide interaction between residues i and j (|j|<=i+1) in antiparallel peptides:
M4=[0, 0, 0, 0, 0, 1, 1;
    0, 0, 0, 0, 1, 1, 1;
    0, 0, 0, 1, 1, 1, 0;
    0, 0, 1, 1, 1, 0, 0;
    0, 1, 1, 1, 0, 0, 0;
    1, 1 ,1, 0, 0, 0, 0;
    1, 1, 0, 0, 0, 0, 0];

% a is the average bond length between two neighboring residues
a=3.8;
% pn is number of peptides included in the calculation. By varying this
% value, one can calculate the free energy surface for oligomers in
% different size simutaneously.
pn=6;
% x is an internal variable for the number of calucation steps
x=0; 
% step is change of constraint parameter C in each step
step=0.1;
%The strength of constraints for each residue in single peptide
Cp(1,1)=0;Cp(1,2)=0;Cp(1,3)=0;Cp(1,4)=0;Cp(1,5)=0;Cp(1,6)=0;Cp(1,7)=0;
Cpre1=1;Cpre2=1;Cpre3=1;Cpre4=1;Cpre5=1;Cpre6=1;Cpre7=1;


for x=2:1:80
    tempF=1000;
for c=1:1:7
    c1=Cpre1;c2=Cpre2;c3=Cpre3;c4=Cpre4;c5=Cpre5;c6=Cpre6;c7=Cpre7;
    % all possible directions for residue constraint evolving 
        switch c
        case 1
                c1=c1+step;
        case 2
                c2=c2+step;
        case 3
                c3=c3+step;
        case 4
                c4=c4+step;
        case 5
                c5=c5+step;
        case 6
                c6=c6+step;
        case 7
                c7=c7+step;
        end

       %Constraint matrix for hexamer
       C=[ c1, 0, 0, 0, 0, 0, 0;
           0, c2, 0, 0, 0, 0, 0;
           0, 0, c3, 0, 0, 0, 0;
           0, 0, 0, c4, 0, 0, 0;
           0, 0, 0, 0, c5, 0, 0;
           0, 0 ,0, 0, 0, c6, 0;
           0, 0, 0, 0, 0, 0, c7];
%C6 is the constraints for hexamer
C6=kron(I6,C);
%B is a constant for controlling peptide collapse
B=0.001;
%Xsich is the inverse monomer correlations for single collapsed stiff chain
Xsich=Xsi+B*I7;
%Xsi0 is the matrix of coefficients for the quadratic terms in Hamitonian H_0
Xsi0=Xsich+C;
%G is the monomer correlations for single collapsed stiff peptide
G=inv(Xsi0);
%Xsich6 is the inverse monomer correlations for hexamer
Xsich6=kron(I6,Xsich);
%G6 is the monomer correlations for hexamer
G6=kron(I6,G);

%s is the average position for each residue
s(7*6,3)=0;
for k=1:3
    for i=1:7*pn
        for j=1:7*pn
            s(i,k)=s(i,k)+G6(i,j)*C6(j,j)*rn(j,k);
        end
    end
end

% interaction matrix constitued by intra-peptide and inter-peptide interaction
% epsilon(i,j)=1 if two residues are in contact; epsilon(i,j)=0 if else.
epsilon=kron(I6,M1)+kron(I6p,M4)+kron(I6a,M3);

% strength parameters for Gaussian-form interaction between residues:
alpha1=3;gamma1=25;
alpha2=0.54;gamma2=9;
alpha3=0.27;gamma3=-6;

%Energy calculations
E=0;
for i=1:7*pn
    for j=i+1:7*pn
        E=E+epsilon(i,j)*WJP(Seq(i),Seq(j))*gamma1/(1+alpha1*(G6(i,i)+G6(j,j)-2*G6(i,j)))^(3/2)*exp(-3/2/a/a*alpha1*(s(i,:)-s(j,:))*(s(i,:)'-s(j,:)')/(1+alpha1*(G6(i,i)+G6(j,j)-2*G6(i,j))));
        E=E+epsilon(i,j)*WJP(Seq(i),Seq(j))*gamma2/(1+alpha2*(G6(i,i)+G6(j,j)-2*G6(i,j)))^(3/2)*exp(-3/2/a/a*alpha2*(s(i,:)-s(j,:))*(s(i,:)'-s(j,:)')/(1+alpha2*(G6(i,i)+G6(j,j)-2*G6(i,j))));
        E=E+epsilon(i,j)*WJP(Seq(i),Seq(j))*gamma3/(1+alpha3*(G6(i,i)+G6(j,j)-2*G6(i,j)))^(3/2)*exp(-3/2/a/a*alpha3*(s(i,:)-s(j,:))*(s(i,:)'-s(j,:)')/(1+alpha3*(G6(i,i)+G6(j,j)-2*G6(i,j))));
    end
end

%Entropy calculations:
S=0;
S=S+3/2*log(det(G))*pn;
for i=1:7*pn
    S=S+3/2*C6(i,i)*G6(i,i);
end
for i=1:7*pn
    for j=1:7*pn
        S=S-3/2/a/a*s(i,:)*Xsich6(i,j)*s(j,:)';
    end
end

%Direction choice of free energy fast descent
if (tempF>E-S)
    tempF=E-S;
    tempc1=c1;tempc2=c2;tempc3=c3;tempc4=c4;tempc5=c5;tempc6=c6;tempc7=c7;
end
    
end

%Free energy calculations
Fp(x)=tempF;
Cp(x,1)=tempc1;Cp(x,2)=tempc2;Cp(x,3)=tempc3;Cp(x,4)=tempc4;Cp(x,5)=tempc5;Cp(x,6)=tempc6;Cp(x,7)=tempc7;
Cpre1=tempc1;Cpre2=tempc2;Cpre3=tempc3;Cpre4=tempc4;Cpre5=tempc5;Cpre6=tempc6;Cpre7=tempc7;
end