Tutorial 4. Fitting of NMR titration series
by Evgenii Kovrigin, 05/17/2011
In this tutorial we will explore fitting of the series of NMR line shapes with different models. We will synthesize data using a complex multi-state model and then fit data with a simpler one as example (we do not doubt the correct model will fit data---what about the wrong one?).
Contents
Clean up workspace (to avoid unwanted effects of old variable assignments)
close all clear all
CHOOSE PARAMETERS OF THE TITRATION SERIES
NOTE: L/R cannot be ZERO. Set some very small number instead!
Rtotal = 1e-3; LRratio_array = [1e-6 0.2 0.4 0.6 0.8 1.0 1.2]; Y_RMSD=1e-7; % ideal line shape Y_offset_percent=0; % to offset figures vertically for easier viewing
CREATE A SERIES OF DATASETS FOR SIMULATION OF LINE SHAPES IN TITRATION
We will use a complex model U-R: Transition A: R+L<=>RL Transition B: R<=>R*
session_name='Tutorial_4';
session=TotalFit(session_name);
Create a simulation series
series_name='test_series'; % arbitrary dataset_class_name='NMRLineShapes1D'; % class name model_name='U_R-model'; model_numeric_solver='analytical'; model_numeric_options='no options'; number_of_datasets=length(LRratio_array); number_of_points_to_simulate=30; X_min=-200; X_max=800; X_column=linspace(X_min, X_max, number_of_points_to_simulate)'; % use ' to transpose to a column-vector session.create_simulation_series(series_name, dataset_class_name, model_name, ... model_numeric_solver, model_numeric_options,... number_of_datasets, X_column); session.show_series_datasets(series_name) % show only datasets of the series session.initialize_relation_matrix(sprintf('%s.all_datasets.txt', session_name)) session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name))
ans = =========== Series "test_series" =============== ----------------------------- Dataset 1 --------------- Priority = 1 Name: test_series-1 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 2 --------------- Priority = 1 Name: test_series-2 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 3 --------------- Priority = 1 Name: test_series-3 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 4 --------------- Priority = 1 Name: test_series-4 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 5 --------------- Priority = 1 Name: test_series-5 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 6 --------------- Priority = 1 Name: test_series-6 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset 7 --------------- Priority = 1 Name: test_series-7 Type: 1D NMR Line shapes Active model: Model 6: "U_R-model" Model description "1D NMR line shape for the U-R model, no temperature/field dependence, analytical" Model handle: line_shape_equations_1D.U_R_model_1D Current solver: analytical Model parameters: 1: Rtotal 2: LRratio 3: log10(K_A) 4: log10(K_B) 5: k_2_A 6: k_2_B 7: w0_R 8: w0_Rstar 9: w0_RL 10: FWHH_R 11: FWHH_Rstar 12: FWHH_RL 13: ScaleFactor Dataset listing is written to "Tutorial_4.all_datasets.txt". Fitting environment saved in Tutorial_4.fitting_environment.txt. ans = New parameter relations (indexed as in all_parameters vector): [1] 1:Rtotal [2] 1:LRratio [3] 1:log10(K_A) [4] 1:log10(K_B) [5] 1:k_2_A [6] 1:k_2_B [7] 1:w0_R [8] 1:w0_Rstar [9] 1:w0_RL [10] 1:FWHH_R [11] 1:FWHH_Rstar [12] 1:FWHH_RL [13] 1:ScaleFactor [14] 2:Rtotal [15] 2:LRratio [16] 2:log10(K_A) [17] 2:log10(K_B) [18] 2:k_2_A [19] 2:k_2_B [20] 2:w0_R [21] 2:w0_Rstar [22] 2:w0_RL [23] 2:FWHH_R [24] 2:FWHH_Rstar [25] 2:FWHH_RL [26] 2:ScaleFactor [27] 3:Rtotal [28] 3:LRratio [29] 3:log10(K_A) [30] 3:log10(K_B) [31] 3:k_2_A [32] 3:k_2_B [33] 3:w0_R [34] 3:w0_Rstar [35] 3:w0_RL [36] 3:FWHH_R [37] 3:FWHH_Rstar [38] 3:FWHH_RL [39] 3:ScaleFactor [40] 4:Rtotal [41] 4:LRratio [42] 4:log10(K_A) [43] 4:log10(K_B) [44] 4:k_2_A [45] 4:k_2_B [46] 4:w0_R [47] 4:w0_Rstar [48] 4:w0_RL [49] 4:FWHH_R [50] 4:FWHH_Rstar [51] 4:FWHH_RL [52] 4:ScaleFactor [53] 5:Rtotal [54] 5:LRratio [55] 5:log10(K_A) [56] 5:log10(K_B) [57] 5:k_2_A [58] 5:k_2_B [59] 5:w0_R [60] 5:w0_Rstar [61] 5:w0_RL [62] 5:FWHH_R [63] 5:FWHH_Rstar [64] 5:FWHH_RL [65] 5:ScaleFactor [66] 6:Rtotal [67] 6:LRratio [68] 6:log10(K_A) [69] 6:log10(K_B) [70] 6:k_2_A [71] 6:k_2_B [72] 6:w0_R [73] 6:w0_Rstar [74] 6:w0_RL [75] 6:FWHH_R [76] 6:FWHH_Rstar [77] 6:FWHH_RL [78] 6:ScaleFactor [79] 7:Rtotal [80] 7:LRratio [81] 7:log10(K_A) [82] 7:log10(K_B) [83] 7:k_2_A [84] 7:k_2_B [85] 7:w0_R [86] 7:w0_Rstar [87] 7:w0_RL [88] 7:FWHH_R [89] 7:FWHH_Rstar [90] 7:FWHH_RL [91] 7:ScaleFactor Summary of datasets with lookup vectors (position of each parameter in all_parameter vector) -------- Dataset 1/test_series-1, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 1 2 3 4 5 6 7 8 9 10 11 12 13 Dataset 2/test_series-2, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 14 15 16 17 18 19 20 21 22 23 24 25 26 Dataset 3/test_series-3, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 27 28 29 30 31 32 33 34 35 36 37 38 39 Dataset 4/test_series-4, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 40 41 42 43 44 45 46 47 48 49 50 51 52 Dataset 5/test_series-5, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 53 54 55 56 57 58 59 60 61 62 63 64 65 Dataset 6/test_series-6, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 66 67 68 69 70 71 72 73 74 75 76 77 78 Dataset 7/test_series-7, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 79 80 81 82 83 84 85 86 87 88 89 90 91
CREATE PARAMETER LINKAGES
parameter_number_array=[ 3 4 5 6 7 8 9 10 11 12 13 ]; session.batch_link_parameters_in_series(series_name, parameter_number_array);
Propagate information about links across all datasets
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name))
Fitting environment saved in Tutorial_4.fitting_environment.txt. ans = New parameter relations (indexed as in all_parameters vector): [1] 1:Rtotal [2] 1:LRratio [3] 1:log10(K_A) | 2:log10(K_A) | 3:log10(K_A) | 4:log10(K_A) | 5:log10(K_A) | 6:log10(K_A) | 7:log10(K_A) [4] 1:log10(K_B) | 2:log10(K_B) | 3:log10(K_B) | 4:log10(K_B) | 5:log10(K_B) | 6:log10(K_B) | 7:log10(K_B) [5] 1:k_2_A | 2:k_2_A | 3:k_2_A | 4:k_2_A | 5:k_2_A | 6:k_2_A | 7:k_2_A [6] 1:k_2_B | 2:k_2_B | 3:k_2_B | 4:k_2_B | 5:k_2_B | 6:k_2_B | 7:k_2_B [7] 1:w0_R | 2:w0_R | 3:w0_R | 4:w0_R | 5:w0_R | 6:w0_R | 7:w0_R [8] 1:w0_Rstar | 2:w0_Rstar | 3:w0_Rstar | 4:w0_Rstar | 5:w0_Rstar | 6:w0_Rstar | 7:w0_Rstar [9] 1:w0_RL | 2:w0_RL | 3:w0_RL | 4:w0_RL | 5:w0_RL | 6:w0_RL | 7:w0_RL [10] 1:FWHH_R | 2:FWHH_R | 3:FWHH_R | 4:FWHH_R | 5:FWHH_R | 6:FWHH_R | 7:FWHH_R [11] 1:FWHH_Rstar | 2:FWHH_Rstar | 3:FWHH_Rstar | 4:FWHH_Rstar | 5:FWHH_Rstar | 6:FWHH_Rstar | 7:FWHH_Rstar [12] 1:FWHH_RL | 2:FWHH_RL | 3:FWHH_RL | 4:FWHH_RL | 5:FWHH_RL | 6:FWHH_RL | 7:FWHH_RL [13] 1:ScaleFactor | 2:ScaleFactor | 3:ScaleFactor | 4:ScaleFactor | 5:ScaleFactor | 6:ScaleFactor | 7:ScaleFactor [14] 2:Rtotal [15] 2:LRratio [16] n/a [17] n/a [18] n/a [19] n/a [20] n/a [21] n/a [22] n/a [23] n/a [24] n/a [25] n/a [26] n/a [27] 3:Rtotal [28] 3:LRratio [29] n/a [30] n/a [31] n/a [32] n/a [33] n/a [34] n/a [35] n/a [36] n/a [37] n/a [38] n/a [39] n/a [40] 4:Rtotal [41] 4:LRratio [42] n/a [43] n/a [44] n/a [45] n/a [46] n/a [47] n/a [48] n/a [49] n/a [50] n/a [51] n/a [52] n/a [53] 5:Rtotal [54] 5:LRratio [55] n/a [56] n/a [57] n/a [58] n/a [59] n/a [60] n/a [61] n/a [62] n/a [63] n/a [64] n/a [65] n/a [66] 6:Rtotal [67] 6:LRratio [68] n/a [69] n/a [70] n/a [71] n/a [72] n/a [73] n/a [74] n/a [75] n/a [76] n/a [77] n/a [78] n/a [79] 7:Rtotal [80] 7:LRratio [81] n/a [82] n/a [83] n/a [84] n/a [85] n/a [86] n/a [87] n/a [88] n/a [89] n/a [90] n/a [91] n/a Summary of datasets with lookup vectors (position of each parameter in all_parameter vector) -------- Dataset 1/test_series-1, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 1 2 3 4 5 6 7 8 9 10 11 12 13 Dataset 2/test_series-2, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 14 15 3 4 5 6 7 8 9 10 11 12 13 Dataset 3/test_series-3, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 27 28 3 4 5 6 7 8 9 10 11 12 13 Dataset 4/test_series-4, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 40 41 3 4 5 6 7 8 9 10 11 12 13 Dataset 5/test_series-5, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 53 54 3 4 5 6 7 8 9 10 11 12 13 Dataset 6/test_series-6, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 66 67 3 4 5 6 7 8 9 10 11 12 13 Dataset 7/test_series-7, model "U_R-model" Parameters : Rtotal LRratio log10(K_A) log10(K_B) k_2_A k_2_B w0_R w0_Rstar w0_RL FWHH_R FWHH_Rstar FWHH_RL ScaleFactor Lookup vector : 79 80 3 4 5 6 7 8 9 10 11 12 13
ASSIGN PARAMETERS
Choose some reasonable parameters (taken from Bill Hawse's spectra of MHC residue Q72)
setup.log10_K_A = log10(1e7) ; % Assume tight binding setup.log10_K_B = log10(3) ; % Assume significant population shift towards non-binding species setup.k_2_A = 10 ; % slow off rate setup.k_2_B = 500 ; % fast isomerization setup.w0_R = 0 ; % resonance frequencies in units 1/s (Hz*2pi) setup.w0_Rstar=100; setup.w0_RL = 500; setup.FWHH_R = 150; % /s (Hz*2pi) setup.FWHH_Rstar = 150 ; setup.FWHH_RL = 190 ; % complex formation increases molecular weight setup.ScaleFactor = 1; % (no meaning in simulation context)
prepare vectors for easier handling:
FAKE_VALUE=1E-6; % individual (unlinked) parameters: to be assigned later ideal_parameters=[ FAKE_VALUE FAKE_VALUE setup.log10_K_A setup.log10_K_B ... setup.k_2_A setup.k_2_B ... setup.w0_R setup.w0_Rstar setup.w0_RL... setup.FWHH_R setup.FWHH_Rstar setup.FWHH_RL... setup.ScaleFactor ] fake_limits=zeros(1,length(ideal_parameters));
ideal_parameters = Columns 1 through 5 0.0000 0.0000 7.0000 0.4771 10.0000 Columns 6 through 10 500.0000 0 100.0000 500.0000 150.0000 Columns 11 through 13 150.0000 190.0000 1.0000
Assign linked parameters (I keep names of function arguments for information)
first_dataset_index=session.dataset_index(series_name,1);
Monte_Carlo_range_min=fake_limits;
Monte_Carlo_range_max=fake_limits;
parameter_limits_min=fake_limits;
parameter_limits_max=fake_limits;
session.assign_parameter_values(first_dataset_index, ideal_parameters, ...
Monte_Carlo_range_min, Monte_Carlo_range_max, parameter_limits_min, parameter_limits_max);
Assign unlinked parameters
Rtotal_vector=ones(1,length(LRratio_array))*Rtotal; % this one is constant for all datasets
Create vectors of ranges for formal assignment (NOTE: they are different length from the above!)
unlinked_fake_ranges=zeros(1,number_of_datasets); MonteCarlo_range_min_array=unlinked_fake_ranges; MonteCarlo_range_max_array=unlinked_fake_ranges; parameter_limits_min_array=unlinked_fake_ranges; parameter_limits_max_array=unlinked_fake_ranges;
Assign Rtotal
parameter_number=1; values_array=Rtotal_vector; session.assign_unlinked_parameter_in_series(... series_name, parameter_number, values_array, ... MonteCarlo_range_min_array, MonteCarlo_range_max_array, ... parameter_limits_min_array, parameter_limits_max_array);
Assign LRratio
parameter_number=2; values_array=LRratio_array; session.assign_unlinked_parameter_in_series(... series_name, parameter_number, values_array, ... MonteCarlo_range_min_array, MonteCarlo_range_max_array, ... parameter_limits_min_array, parameter_limits_max_array);
CALCULATE MODEL VALUES
Run a simulation
X_standard_dev=0; % no need to have errors of Y
session.simulate_series(series_name, X_standard_dev, Y_RMSD);
Plot the simulation result
figure_handle1=session.plot_1D_series(series_name, 'Simulation 4', Y_offset_percent);
output simulation
figure_handle=figure_handle1; folder_name=session_name; figure_name=series_name; results_output.output_figure(figure_handle, folder_name, figure_name)
PREPARE FOR FITTING WITH SIMPLE U MODEL
Save original simulation session
session_simulation=session.make_a_copy('Tutorial_4_simulation');
Set new session name
session_name='Tutorial_4_U_fit';
session.set_name(session_name);
Switch model
new_model_name='U-model'; model_numeric_solver='analytical'; model_numeric_options='no options'; session.change_model_in_series(series_name, new_model_name,... model_numeric_solver, model_numeric_options) session.initialize_relation_matrix(sprintf('%s.all_datasets.txt', session_name)) session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name))
Dataset "test_series-1" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-2" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-3" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-4" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-5" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-6" has the model switched to "U-model". All parameter links and values are reset. Dataset "test_series-7" has the model switched to "U-model". All parameter links and values are reset. Dataset listing is written to "Tutorial_4_U_fit.all_datasets.txt". Fitting environment saved in Tutorial_4_U_fit.fitting_environment.txt. ans = New parameter relations (indexed as in all_parameters vector): [1] 1:Rtotal [2] 1:LRratio [3] 1:log10(K_A) [4] 1:k_2_A [5] 1:w0_1 [6] 1:w0_2 [7] 1:FWHH_1 [8] 1:FWHH_2 [9] 1:ScaleFactor [10] 2:Rtotal [11] 2:LRratio [12] 2:log10(K_A) [13] 2:k_2_A [14] 2:w0_1 [15] 2:w0_2 [16] 2:FWHH_1 [17] 2:FWHH_2 [18] 2:ScaleFactor [19] 3:Rtotal [20] 3:LRratio [21] 3:log10(K_A) [22] 3:k_2_A [23] 3:w0_1 [24] 3:w0_2 [25] 3:FWHH_1 [26] 3:FWHH_2 [27] 3:ScaleFactor [28] 4:Rtotal [29] 4:LRratio [30] 4:log10(K_A) [31] 4:k_2_A [32] 4:w0_1 [33] 4:w0_2 [34] 4:FWHH_1 [35] 4:FWHH_2 [36] 4:ScaleFactor [37] 5:Rtotal [38] 5:LRratio [39] 5:log10(K_A) [40] 5:k_2_A [41] 5:w0_1 [42] 5:w0_2 [43] 5:FWHH_1 [44] 5:FWHH_2 [45] 5:ScaleFactor [46] 6:Rtotal [47] 6:LRratio [48] 6:log10(K_A) [49] 6:k_2_A [50] 6:w0_1 [51] 6:w0_2 [52] 6:FWHH_1 [53] 6:FWHH_2 [54] 6:ScaleFactor [55] 7:Rtotal [56] 7:LRratio [57] 7:log10(K_A) [58] 7:k_2_A [59] 7:w0_1 [60] 7:w0_2 [61] 7:FWHH_1 [62] 7:FWHH_2 [63] 7:ScaleFactor Summary of datasets with lookup vectors (position of each parameter in all_parameter vector) -------- Dataset 1/test_series-1, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 1 2 3 4 5 6 7 8 9 Dataset 2/test_series-2, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 10 11 12 13 14 15 16 17 18 Dataset 3/test_series-3, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 19 20 21 22 23 24 25 26 27 Dataset 4/test_series-4, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 28 29 30 31 32 33 34 35 36 Dataset 5/test_series-5, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 37 38 39 40 41 42 43 44 45 Dataset 6/test_series-6, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 46 47 48 49 50 51 52 53 54 Dataset 7/test_series-7, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 55 56 57 58 59 60 61 62 63
CREATE PARAMETER LINKAGES
parameter_number_array=[ 3 4 5 6 7 8 9 ]; session.batch_link_parameters_in_series(series_name, parameter_number_array);
Propagate information about links across all datasets
session.generate_fitting_environment(sprintf('%s.fitting_environment.txt', session_name))
Fitting environment saved in Tutorial_4_U_fit.fitting_environment.txt. ans = New parameter relations (indexed as in all_parameters vector): [1] 1:Rtotal [2] 1:LRratio [3] 1:log10(K_A) | 2:log10(K_A) | 3:log10(K_A) | 4:log10(K_A) | 5:log10(K_A) | 6:log10(K_A) | 7:log10(K_A) [4] 1:k_2_A | 2:k_2_A | 3:k_2_A | 4:k_2_A | 5:k_2_A | 6:k_2_A | 7:k_2_A [5] 1:w0_1 | 2:w0_1 | 3:w0_1 | 4:w0_1 | 5:w0_1 | 6:w0_1 | 7:w0_1 [6] 1:w0_2 | 2:w0_2 | 3:w0_2 | 4:w0_2 | 5:w0_2 | 6:w0_2 | 7:w0_2 [7] 1:FWHH_1 | 2:FWHH_1 | 3:FWHH_1 | 4:FWHH_1 | 5:FWHH_1 | 6:FWHH_1 | 7:FWHH_1 [8] 1:FWHH_2 | 2:FWHH_2 | 3:FWHH_2 | 4:FWHH_2 | 5:FWHH_2 | 6:FWHH_2 | 7:FWHH_2 [9] 1:ScaleFactor | 2:ScaleFactor | 3:ScaleFactor | 4:ScaleFactor | 5:ScaleFactor | 6:ScaleFactor | 7:ScaleFactor [10] 2:Rtotal [11] 2:LRratio [12] n/a [13] n/a [14] n/a [15] n/a [16] n/a [17] n/a [18] n/a [19] 3:Rtotal [20] 3:LRratio [21] n/a [22] n/a [23] n/a [24] n/a [25] n/a [26] n/a [27] n/a [28] 4:Rtotal [29] 4:LRratio [30] n/a [31] n/a [32] n/a [33] n/a [34] n/a [35] n/a [36] n/a [37] 5:Rtotal [38] 5:LRratio [39] n/a [40] n/a [41] n/a [42] n/a [43] n/a [44] n/a [45] n/a [46] 6:Rtotal [47] 6:LRratio [48] n/a [49] n/a [50] n/a [51] n/a [52] n/a [53] n/a [54] n/a [55] 7:Rtotal [56] 7:LRratio [57] n/a [58] n/a [59] n/a [60] n/a [61] n/a [62] n/a [63] n/a Summary of datasets with lookup vectors (position of each parameter in all_parameter vector) -------- Dataset 1/test_series-1, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 1 2 3 4 5 6 7 8 9 Dataset 2/test_series-2, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 10 11 3 4 5 6 7 8 9 Dataset 3/test_series-3, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 19 20 3 4 5 6 7 8 9 Dataset 4/test_series-4, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 28 29 3 4 5 6 7 8 9 Dataset 5/test_series-5, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 37 38 3 4 5 6 7 8 9 Dataset 6/test_series-6, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 46 47 3 4 5 6 7 8 9 Dataset 7/test_series-7, model "U-model" Parameters : Rtotal LRratio log10(K_A) k_2_A w0_1 w0_2 FWHH_1 FWHH_2 ScaleFactor Lookup vector : 55 56 3 4 5 6 7 8 9
ASSIGN STARTING PARAMETERS
Binding constant
Log10_K_A=log10(1e6); Log10_K_A_AbsMin=log10(1e1); Log10_K_A_AbsMax=log10(1e10);
Dissociation rate constant
k_2_A=100; k_2_A_AbsMin=0.1; k_2_A_AbsMax=3000;
Frequency of R
w0_R=50;
w0_R_AbsMin=X_min; % set these limits to the spectral window
w0_R_AbsMax=X_max;
Frequency of RL
w0_RL=500;
w0_RL_AbsMin=X_min; % set these limits to the spectral window
w0_RL_AbsMax=X_max;
Line width of R
FWHH_R=100; FWHH_R_AbsMin=10; FWHH_R_AbsMax=500;
Line width of RL
FWHH_RL=100; FWHH_RL_AbsMin=10; FWHH_RL_AbsMax=500;
ScaleFactor
ScaleFactor=1; ScaleFactor_AbsMin=0.1; ScaleFactor_AbsMax=10;
Create a vector of parameter values
starting_parameters=[ FAKE_VALUE FAKE_VALUE ... Log10_K_A k_2_A ... w0_R w0_RL FWHH_R FWHH_RL ScaleFactor ];
Create vectors of absolute limits. Parameters, which will not be optimized do not need meaningful limits (use 0).
min_limits=[ 0 0 Log10_K_A_AbsMin k_2_A_AbsMin w0_R_AbsMin w0_RL_AbsMin FWHH_R_AbsMin FWHH_RL_AbsMin ScaleFactor_AbsMin ]; max_limits=[ 0 0 Log10_K_A_AbsMax k_2_A_AbsMax w0_R_AbsMax w0_RL_AbsMax FWHH_R_AbsMax FWHH_RL_AbsMax ScaleFactor_AbsMax ];
Assign linked parameters (I keep names of function arguments for information) I am formally setting Monte-Carlo range to the same arrays as absolute limits. In reality, they will not be used because we are using the default mode vary_starting_parameters='no' meaning that starting parameters are not varied for new Monte-Carlo runs.
first_dataset_index=session.dataset_index(series_name,1);
Monte_Carlo_range_min=min_limits;
Monte_Carlo_range_max=max_limits;
parameter_limits_min=min_limits;
parameter_limits_max=max_limits;
session.assign_parameter_values(first_dataset_index, starting_parameters, ...
Monte_Carlo_range_min, Monte_Carlo_range_max, parameter_limits_min, parameter_limits_max);
Assign unlinked parameters Create vectors of ranges for formal assignment (different length from the above!)
unlinked_fake_ranges=zeros(1,number_of_datasets);
Assign Rtotal --- this will not be fit so no need for ranges
parameter_number=1; values_array=Rtotal_vector; % this remains the same MonteCarlo_range_min_array=unlinked_fake_ranges; MonteCarlo_range_max_array=unlinked_fake_ranges; parameter_limits_min_array=unlinked_fake_ranges; parameter_limits_max_array=unlinked_fake_ranges; session.assign_unlinked_parameter_in_series(... series_name, parameter_number, values_array, ... MonteCarlo_range_min_array, MonteCarlo_range_max_array, ... parameter_limits_min_array, parameter_limits_max_array);
Assign LRratio
parameter_number=2; values_array=LRratio_array; session.assign_unlinked_parameter_in_series(... series_name, parameter_number, values_array, ... MonteCarlo_range_min_array, MonteCarlo_range_max_array, ... parameter_limits_min_array, parameter_limits_max_array);
Inspect initial guess
figure_handle1=session.plot_1D_series(series_name, 'U_fitting', Y_offset_percent);
Good terrible guess.
Save initial state of fitting session on a disk for future inspection and reuse
session.export_results(sprintf('%s.initial_guess', session_name), {figure_handle1})
Output current session state, statistics and graphs into "Tutorial_4_U_fit.initial_guess/". Adjustable parameters and cumulative dataset statistics is written into "Tutorial_4_U_fit.initial_guess/adjustable_parameters_only.txt". Statistics on per-dataset basis is written into "Tutorial_4_U_fit.initial_guess/all_datasets.txt". Writing figures... Figure 1: small PNG image, large PNG image, EPS image, MATLAB figure file - written Done
FIT
Set which parameters of the model will be fixed and variable
parameter_modes=[ 0 0 1 1 1 1 1 1 0]; session.set_parameters_modes(first_dataset_index, parameter_modes);
Select fitting algorithm. Choices: 'simplex', 'MultiStart' 'GlobalSearch', 'DirectSearch'. Do not forget to start your parallel processors: matlabpool open
fit_only_algorithm='GlobalSearch';
Report on start-up and show variable parameter names
fprintf(1,'\n\n==============\nFitting Algorithm: %s\nVariable parameters:', fit_only_algorithm);
session.fit_parameter_names{1:end}
============== Fitting Algorithm: GlobalSearch Variable parameters: ans = 1:log10(K_A) | 2:log10(K_A) | 3:log10(K_A) | 4:log10(K_A) | 5:log10(K_A) | 6:log10(K_A) | 7:log10(K_A) ans = 1:k_2_A | 2:k_2_A | 3:k_2_A | 4:k_2_A | 5:k_2_A | 6:k_2_A | 7:k_2_A ans = 1:w0_1 | 2:w0_1 | 3:w0_1 | 4:w0_1 | 5:w0_1 | 6:w0_1 | 7:w0_1 ans = 1:w0_2 | 2:w0_2 | 3:w0_2 | 4:w0_2 | 5:w0_2 | 6:w0_2 | 7:w0_2 ans = 1:FWHH_1 | 2:FWHH_1 | 3:FWHH_1 | 4:FWHH_1 | 5:FWHH_1 | 6:FWHH_1 | 7:FWHH_1 ans = 1:FWHH_2 | 2:FWHH_2 | 3:FWHH_2 | 4:FWHH_2 | 5:FWHH_2 | 6:FWHH_2 | 7:FWHH_2
Fit without error determination
tic session.fit_only(fit_only_algorithm); toc
GlobalSearch stopped because it analyzed all the trial points. All 63 local solver runs converged with a positive local solver exit flag. Elapsed time is 70.166584 seconds.
Report on finishing
session.show_fit_statistics() session.transfer_all_parameters_to_individual_datasets();
ans = Fitting results for "Tutorial_4_U_fit" ========================== Fitting method: GlobalSearch Number of active datasets (non-zero priority): 7 Reporting only on variable parameters (dataset-#:parameter-name): (# in the fit-vector) [min value, max value] for 95% confidence interval -------------------------- (1) name: " 1:log10(K_A) | 2:log10(K_A) | 3:log10(K_A) | 4:log10(K_A) | 5:log10(K_A) | 6:log10(K_A) | 7:log10(K_A)" (1) value= 4.345454e+00 [ 0.000000e+00 0.000000e+00 ] (2) name: " 1:k_2_A | 2:k_2_A | 3:k_2_A | 4:k_2_A | 5:k_2_A | 6:k_2_A | 7:k_2_A" (2) value= 4.269107e+01 [ 0.000000e+00 0.000000e+00 ] (3) name: " 1:w0_1 | 2:w0_1 | 3:w0_1 | 4:w0_1 | 5:w0_1 | 6:w0_1 | 7:w0_1" (3) value= 7.525018e+01 [ 0.000000e+00 0.000000e+00 ] (4) name: " 1:w0_2 | 2:w0_2 | 3:w0_2 | 4:w0_2 | 5:w0_2 | 6:w0_2 | 7:w0_2" (4) value= 5.281800e+02 [ 0.000000e+00 0.000000e+00 ] (5) name: " 1:FWHH_1 | 2:FWHH_1 | 3:FWHH_1 | 4:FWHH_1 | 5:FWHH_1 | 6:FWHH_1 | 7:FWHH_1" (5) value= 1.358092e+02 [ 0.000000e+00 0.000000e+00 ] (6) name: " 1:FWHH_2 | 2:FWHH_2 | 3:FWHH_2 | 4:FWHH_2 | 5:FWHH_2 | 6:FWHH_2 | 7:FWHH_2" (6) value= 1.131700e+02 [ 0.000000e+00 0.000000e+00 ] Total sum of squares (SS) = 7.4495e-06 Statistics per dataset: ------------ (1) "test_series-1": Priority=1, R^2=0.991, 15% of total SS, Model= "U-model" (2) "test_series-2": Priority=1, R^2=0.994, 5.1% of total SS, Model= "U-model" (3) "test_series-3": Priority=1, R^2=0.968, 14% of total SS, Model= "U-model" (4) "test_series-4": Priority=1, R^2=0.946, 18% of total SS, Model= "U-model" (5) "test_series-5": Priority=1, R^2=0.983, 9.3% of total SS, Model= "U-model" (6) "test_series-6": Priority=1, R^2=0.980, 23% of total SS, Model= "U-model" (7) "test_series-7": Priority=1, R^2=0.987, 15% of total SS, Model= "U-model"
Inspect final result
figure_handle2=session.plot_1D_series(series_name, 'U_final', Y_offset_percent);
Export final result
session.export_results(sprintf('%s.Fit_Only.%s', session_name, fit_only_algorithm), {figure_handle2})
Output current session state, statistics and graphs into "Tutorial_4_U_fit.Fit_Only.GlobalSearch/". Adjustable parameters and cumulative dataset statistics is written into "Tutorial_4_U_fit.Fit_Only.GlobalSearch/adjustable_parameters_only.txt". Statistics on per-dataset basis is written into "Tutorial_4_U_fit.Fit_Only.GlobalSearch/all_datasets.txt". Writing figures... Figure 1: small PNG image, large PNG image, EPS image, MATLAB figure file - written Done
CONCLUSION
In slow exchange, two-state model easily can be fit to the three-state data: we are not sensitive to the model structure enough.