Fitting Inversion Magnetization Transfer NMR Experiments

DESCRIPTION

	This notebook will simultaneously fit the time-dependent magnetization transfers of Pro 1H or 13C cis (Mct) and trans (Mtt) NMR signals as a function of the inversion transfer time (t) in Inversion Magnetization Transfer experiments.  In the following pulse sequence, the 13C or 1H cis resonance is selectively inverted using a shaped pulse.  Its recovery during t is determined by its intrinsic T1c, magnetization transfer to and from the trans resonance, and the T1t of the trans resonance.  Similarly, the steady-state magnetization of the trans resonance is modulated by exchange with the cis resonance:
	
	               p(x)sel ------------t-----------p/2(x,y.-x.-y)---acquire
	               
The resonance of the trans isomer of Pro (Mtt) shows the following time dependence:
	               
         Mtt = (c1)(tt)(l1+1/t1c)exp(l1*t) + (c2)(t2)(l2+1/t1c)exp(l2t) + Mtinf
	 
The resonance of the cis isomer of Pro (Mtc) shows the following time dependence:
	 
	                            Mtc = (c1)exp(l1*t) + (c2)exp(l2t) + Mcinf
	                            
T1c and T1t are the longitudinal relaxation times the resonances would have in the absence of exchange.

tc and tt are the lifetimes of the cis and trans conformers and kct and ktc are the corresponding rate constants.

t1c and t1t are the effective relaxation times of the cis and trans resonances when relaxation and exchange are both occuring and are defined below in terms of T1c and tc, T1t and tt.

l1 and l2 are related to the time constants tc, tt, t1c, and t1t, and are defined below.

c1 and c2 are defined below.

The user must enter Mcinf and Mtinf, the magnetization measured after 5 T1 periods for the cis and trans resonances.

The program calculates tt from tc, Mcinf, and Mtinf as: tt = tc * Mtinf/Mcinf.

The notebook will also generate graphs of the data and the fitted curve and  statistical information on the goodness of the fit. A table of points can be produced to permit export of the theoretical curve.	
	                            
References:

Alger and Prestegard (1977) J. Magn. Reson. 27, 137-41.
Mariappan and Rabenstein (1992) J. Magn. Reson. 100, 183-8.

Thanks to Maxim A. Dubinnyi for advice on how to fit multiple data sets.	                        	                                        			        		
    			   					                                        			        		INSTRUCTIONS

	Hit the down arrow to move from cell to cell and press "enter" (not return) to execute each calculation. Follow the example below to see the fitting results for one set of data for hydroxy-Pro.
	
	The user may enter data manually, in the form of {Transfer Time, Intensity} pairs, or read in a data files. For the latter, set the directory containing the data in the first line below.  If you are not reading in files, move the cursor to the 4th line ("list2") and replace the pairs of points in "list2" with your own data.  Do the same with "list3". Alternatively, press enter on the line "list2" and follow the fitting of those data.
		
	In the first 3 lines below we set a directory, read in a data file, and check the number of points in the file.

SetDirectory["~/Documents/Craig/"];

list2=ReadList["HydroxyProtrans.txt",Number,RecordLists->True];
list3=ReadList["HydroxyProcis.txt",Number,RecordLists->True];

totalres2=Length[list2];
totalres3=Length[list3]

Users: Manual data entry begins here:

list2:={{0.001,111.31},{0.05,108.11},{0.2,101.40},{0.4,93.99},{0.8,85.91},{1.3,83.75},{1.9,89.10},{2.5,93.61},{3.0,95.39},{5.0,103.29},{7.0,106.93},{10.0,108.45}}

list3:={{0.001,-29.18},{0.05,-26.84},{0.2,-13.33},{0.4,-6.77},{0.8,10.00},{1.3,20.30},{1.9,24.06},{2.5,27.66},{3.0,28.22},{5.0,34.68},{7.0,33.28},{10.0,35.32}}

lp2=ListPlot[list2,PlotStyle -> {PointSize[0.02]}, FrameStyle->{Thickness[0.01]},PlotRange->{{0,15},{70,120}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time","Intensity"},Axes->False]

lp3=ListPlot[list3,PlotStyle -> {PointSize[0.02]}, FrameStyle->{Thickness[0.01]},PlotRange->{{0,15},{-75,50}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time","Intensity"},Axes->False]

CombinedData = Join[Insert[#,1,2]&/@list2, Insert[#,2,2]&/@list3];
CombinedData//TableForm

<<"NonlinearRegression`"

Clear[T1t,T1c,taut,tauc,M0t,M0c,Mtinf,Mcinf,tau1t,tau1c,lam1,lam2,c1,c2,Mct,Mtt,Keq,kct,ktc,param,NLR]

Users: Enter your values for Mcinf and Mtinf below:

Mcinf:=35.32

Mtinf:=108.45

taut:=((tauc)*(Mtinf))/(Mcinf)

kct:=1/tauc

ktc:=1/taut

Keq:=(Mtinf)/(Mcinf)

tau1c:=((T1c)*(tauc))/((tauc)+(T1c))

tau1t:=((T1t)*(taut))/((taut)+(T1t))

lam1:=(0.5)*(-(((1)/(tau1c))+((1)/(tau1t)))+((((((1)/(tau1c))+((1)/(tau1t)))^(2))-(4)*((((1)/(tau1c))*((1)/(tau1t)))-(((1)/(tauc))*((1)/(taut))))))^(0.5))

lam2:=(0.5)*(-(((1)/(tau1c))+((1)/(tau1t)))-((((((1)/(tau1c))+((1)/(tau1t)))^(2))-(4)*(((((1)/(tau1c)))*(((1)/(tau1t))))-(((1)/(tauc))*((1)/(taut))))))^(0.5))

c2:=((1)/((taut)*((lam1)-(lam2))))*(((taut)*((lam1)+(((1)/(tau1c))))*((M0c)-(Mcinf)))+((Mtinf)-(M0t)))

c1:=M0c-Mcinf-c2

Mtt:=((c1)*(taut)*((lam1)+(((1)/(tau1c))))*(Exp[(lam1)*(x)]))+((c2)*(taut)*((lam2)+(((1)/(tau1c))))*(Exp[(lam2)*(x)]))+Mtinf

Mct:=((c1)*(Exp[(lam1)*(x)]))+((c2)*(Exp[(lam2)*(x)]))+(Mcinf)

CombinedMctMtt :=
  Which[MDL == 1, Evaluate@Mtt,
        MDL == 2, Evaluate@Mct,
        True,     0];

taut

tau1c

tau1t

1/tau1c

1/tau1t

Users: You may want to adjust the best guesses for the parameters, the upper and lower limits, and the number of interations below:

NLR=NonlinearRegress[CombinedData,CombinedMctMtt, {{T1c,3,0.5,10},{T1t,2,0.5,10},{tauc,1,0.1,10},{M0t,111,80,200},{M0c,-29,-50,-20}}, {x,MDL},MaxIterations->10000,ShowProgress ->False ]

param=BestFitParameters/.NLR

T1c=T1c/.param
T1t=T1t/.param
tauc=tauc/.param
M0t=M0t/.param
M0c=M0c/.param

Plot[Mtt, {x,0,15},PlotStyle -> {Thickness[0.01]}, FrameStyle->{Thickness[0.010]},PlotRange->{{0,15},{70,120}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time (sec)","Intensity"},Axes->False]

Show[%,lp2]

Plot[Mct, {x,0,15},PlotStyle -> {Thickness[0.01]}, FrameStyle->{Thickness[0.010]},PlotRange->{{0,15},{-75,50}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time (sec)","Intensity"},Axes->False]

Show[%,lp3]

Following are the rate constants in sec-1 and the equilibrium constant:

kct

ktc

Keq

The next line calculates points from the theoretical curve in case you want to graph the curve in another application:

theorMtt=Table[{x,Mtt},{x,0,15,0.01}]

theorMct=Table[{x,Mct},{x,0,15,0.01}]

The next line exports a file containing the theoretical points just calculated.

Export["OHPro1.txt",theorMtt,"Table"]

Export["OHPro2.txt",theorMct,"Table"]

Just to confirm, here are graphs of the theoretical points and the data.

ListPlot[theorMtt,PlotStyle -> {PointSize[0.01]}, FrameStyle->{Thickness[0.01]},PlotRange->{{0,15},{70,+120}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time (sec)","Intensity"},Axes->False]

Show[%,lp2]

ListPlot[theorMct,PlotStyle -> {PointSize[0.01]}, FrameStyle->{Thickness[0.01]},PlotRange->{{0,15},{-75,+50}}, Frame->True, RotateLabel->False,BaseStyle->{14,FontFamily->"Times"},FrameLabel->{"Transfer Time (sec)","Intensity"},Axes->False]

Show[%,lp3]

Clear[T1t,T1c,taut,tauc,M0t,M0c,Mtinf,Mcinf,tau1t,tau1c,lam1,lam2,c1,c2,Mct,Mtt,Keq,kct,ktc]

END