I recently needed to fit the extensible worm-like chain (EWLC) model to several force-extension curves from some of my single molecule experiments. In these experiments, a DNA/RNA hybrid was stretched using optical tweezers and I wished to estimate the persistence length, the contour length and the stretch modulus.
As the EWLC is an implicit function, this requires implicit fitting. I have implemented this in iGor Pro and will provide the code below so that you don’t have to waste time implementing it yourself.
It is assumed that you have the extension data in a wave called “[prefix]_y_real_dist” and the force data in a wave called “[prefix]_y_force”. If not, update the code to match your wave names.
First the actual implicit function, it takes three arguments:
- w: wave of parameters for the EWLC:
w[0] = kB
w[1] = T
w[2] = Lp
w[3] = Stretch modulus
w[4] = Contour length
w[5] = Extension offset (hold at zero if you know you have zero initial extension) - x: Extension of tether (not a wave, a single value)
- y: Force at extension (not a wave, a single value)
Function FitEWLC(w,x,y): FitFunc
Wave w
Variable x
Variable y
return (w[0]*w[1])/w[2] * ((1/(4*(1-(x+w[5])/w[4] + y/w[3])^2)) - 1/4 + (x+w[5])/w[4] - y/w[3]) - y
End
For performing the actual fitting I have written a small function, DoFitEWLC, which lets you control a few things about the fitting.
DoFitEWLC takes four arguments:
- prefix: I use this because I have many waves with similar names except for the prefix. The prefix is a string.
- forceLimit: If you don’t want to fit your entire dataset, enter the upper limit of the force
- offset: you guess of the initial extension
- holdStr: If you want to hold any of the EWLC parameters, you probably want to, enter the hold string here. Example, if you want to hold the boltzmann constant and temperature the holdStr should be “110000″
Function DoFitEWLC(prefix, forceLimit, offset, holdStr)
//--
String prefix
Variable forceLimit
Variable offset
String holdStr
//--
// Declare wave names with source data
String real_dist
String force
sprintf real_dist, "%s_y_real_dist" prefix
sprintf force, "%s_y_force" prefix
wave w_real_dist = $real_dist
wave w_force = $force
if (!WaveExists(w_real_dist))
DoAlert/T="DataError" 0, "Extension:\rUnable to locate wave for trap1 Y-distance data\rLooked for wave named: \""+real_dist + "\""
return -1
endif
if (!WaveExists(w_force))
DoAlert/T="DataError" 0, "Force:\rUnable to locate wave for trap1 Y-force data\rLooked for wave named: \""+force + "\""
return -1
endif
// Find the range in which we want to fit
Variable n = numpnts(w_force)
Variable i = 0
Variable curr_F = 0
do
curr_F = w_force[i]
i += 1
// If limit is higher than the highest force - use full wave
if (i > n)
i = n
break
endif
while (curr_F < forceLimit)
// Declare waves for fit
String ewlcFitY
String ewlcFitX
sprintf ewlcFitY, "%s_EWLC_Y" prefix
sprintf ewlcFitX, "%s_EWLC_X" prefix
wave w_ewlcFitY = $ewlcFitY
wave w_ewlcFitX = $ewlcFitX
// Declare waves to hold temporary force and disance
String tmp_real_dist = "EWLC_tmp_dist"
String tmp_force = "EWLC_tmp_force"
Wave w_tmp_real_dist = $tmp_real_dist
Wave w_tmp_force = $tmp_force
// Make waves to hold fit values
Duplicate/D/O/R=(0,i) w_real_dist, $ewlcFitX, $ewlcFitY
// Copy the desired range for force and extension
Duplicate/D/O/R=(0,i) w_real_dist, $tmp_real_dist
Duplicate/D/O/R=(0,i) w_force, $tmp_force
// Make wave for parameters
Make/D/O ewlcCoefs={(1.38e-23*1e12*1e9),293,50,1000,300,offset}
FuncFit/ODR=3/H=holdStr FitEWLC, ewlcCoefs /X={$tmp_real_dist,$tmp_force}/XD={$ewlcFitX,$ewlcFitY}
End
The initial guesses for the fitting parameters are:
kB: 1.38e-23 J/K
T: 293 K
Lp: 50 nm
K: 1000 pN
Lc: 300 nm
So, does it work? I simulated data for a tether with a persistence length of 50 nm, a contour length of 300 nm and a stretch modulus of 1000 pN. The force data is in wave “NT_y_force” and extension data is in wave “NT_y_real_dist”. Issuing the following command
DoFitEWLC("NT",30,0,"110001")
Will produce the following output which shows that all parameters were estimated accurately:
Fit converged properly
ewlcCoefs={0.0138,293,50,1000,300,0}
V_chisq= 4.75695e-06;V_npnts= 5987;V_numNaNs= 0;V_numINFs= 0;
W_sigma={0,0,6.07e-05,0.00142,1.4e-05,0}
Coefficient values ± one standard deviation
w_0 =0.0138 ± 0
w_1 =293 ± 0
w_2 =50 ± 6.07e-05
w_3 =1000 ± 0.00142
w_4 =300 ± 1.4e-05
w_5 =0 ± 0
Please note that the holdStr limits the number of free parameters to three and that we only fit at forces up to 30pN.
The figure below shows the simulated data (black) and the fitted data (red).
