Estimate Parameters of a G protein Model

Overview

About the Example Model

This example illustrates parameter estimation using time-course data from one experiment, using the sbioparamestim function. For information on all available parameter estimation and population fitting techniques, see .

This example uses the model described in Model of the Yeast Heterotrimeric G Protein Cycle to illustrate parameter estimation.

This table lists the reactions used to model the G protein cycle and the corresponding rate parameters (rate constants) for each mass action reaction. For reversible reactions, the forward rate parameter is listed first.

No.NameReaction1Rate Parameters
1Receptor-ligand interactionL + R <-> RLkRL, kRLm
2Heterotrimeric G protein formationGd + Gbg -> GkG1
3G protein activationRL + G -> Ga + Gbg + RLkGa
4Receptor synthesis and degradationR <-> nullkRdo, kRs
5Receptor-ligand degradationRL -> nullkRD1
6G protein inactivationGa -> GdkGd
1 Legend of species: L = ligand (alpha factor), R = alpha-factor receptor, Gd = inactive G-alpha-GDP, Gbg = free levels of G-beta:G-gamma complex, G = inactive Gbg:Gd complex, Ga = active G-alpha-GTP

About the Example

The study used to build the example model (Yi et al., 2003) reported the estimated value of parameter kGd as 0.11 for the wild-type strain. In Calculate Sensitivities, the analysis showed that Ga is sensitive to parameters kGd, kRs, kRD1, and kGa.

This example shows:

  • How to estimate the parameter kGd and determine its effect on the model

  • How to estimate parameters kGd, kRs, kRD1, and kGa to obtain a better fit to the experimental data

Loading the Example Model

The gprotein.sbproj project contains a model for the wild-type strain (stored in variable m1). Load the G protein model for the wild-type strain:

sbioloadproject gprotein

The m1 model object appears in the MATLAB® Workspace.

Defining Experimental Data

The study used for this example (Yi et al., 2003) reports, in a plot, the experimental data as the fraction of active G protein. Store the time data for the experimental results in a variable, tExpt, and store the values for the fraction of active G protein in a variable, GaFracExpt:

tExpt  = [0 10 30 60 110 210 300 450 600]';
GaFracExpt = [0 0.35 0.4 0.36 0.39 0.33 0.24 0.17 0.2]';

    Note:   For this simple example, you stored the experimental data in a variable in the MATLAB Workspace by typing the data values. However, typically, you import larger data sets into a MATLAB variable. For more information about importing data into variables, see Methods for Importing Data.

Simulating the G Protein Model

  1. Simulate the model and return the results to a SimData object:

    simDataObj = sbiosimulate(m1);
    
  2. Retrieve the time and state data for the GaFrac parameter:

    [tOrig, GaFracOrig] = selectbyname(simDataObj,'GaFrac');

Calculating R2 for the G Protein Model

R2 is the square of the correlation between the response values and the predicted response values. Therefore, R2 measures how successful the fit is in explaining the variation of the data.

  1. Calculate the sum of squares about the mean (SST):

    sst = norm(GaFracExpt - mean(GaFracExpt))^2;
    
  2. Interpolate the data to get time points that match the time points in the experimental data using the pchip interpolation method:

    GaFracResampled = interp1(tOrig, GaFracOrig, tExpt, 'pchip');
    
  3. Calculate the sum of squares due to error (SSE):

    sse = norm(GaFracExpt - GaFracResampled)^2;
    
  4. Calculate the R2 value for the simulation data before parameter estimation:

    rSquareOrig = 1-sse/sst
    
    rSquareOrig =
    
        0.8968

For more information about the functions used here, see the norm and interp1 reference pages.

Plotting the Experimental Results and Simulation Data

  1. Plot the experimental data for active G protein:

    plot(tExpt, GaFracExpt, 'ro');
    title('Variation of G Protein');
    xlabel('Time (sec)');
    ylabel('Active Fraction of G Protein');
    legend('Experiment');
  2. Plot the simulation data in the same plot:

    hold on;
    plot(tOrig, GaFracOrig);
    legendText = {'Experiment', sprintf('Original R^2=%4.2f',...
                   rSquareOrig)};
    legend(legendText{:});

      Note:   Leave this figure window open so you can use it to plot and compare results of using the estimated parameters later in this example.

Estimating the kGd Parameter in the G Protein Model

The study used to build the G protein model reported an estimated value of 0.11 for the parameter kGd in the wild-type strain (Yi et al., 2003). This example estimates the value of kGd.

  1. Create a variable for the parameter to estimate. Also create a variable for the parameter corresponding to the experimental data to which you are fitting:

    paramToEst = sbioselect(m1, 'Name', 'kGd');
    GaFrac = sbioselect(m1, 'Name', 'GaFrac');
    
  2. Specify plotting of each iteration of the parameter estimation to see how optimization is progressing:

    opt = optimset('PlotFcns',@optimplotfval,'MaxIter',15);
    
  3. Use the current value of the kGd parameter in the model as the starting value for optimization:

    [estValues1, result1] = sbioparamestim(m1, tExpt, GaFracExpt, ...
                             GaFrac, paramToEst, {}, {'fminsearch',opt});

      Note:   Close this figure before proceeding with the example.

Simulating and Plotting Results Using the Estimated Parameter

Use the estimated value of the kGd parameter to see how it affects simulation results.

  1. Use a variant to store the estimated value of kGd:

    estvarObj = addvariant (m1, 'Optimized kGd');
    addcontent(estvarObj, {'parameter', 'kGd', 'Value', estValues1});
  2. Apply the value stored in the variant, simulate the model, and return the results:

    simDataObj1 = sbiosimulate(m1, estvarObj );
    [t1, GaFrac1] = selectbyname(simDataObj1,'GaFrac');
  3. Calculate the R2 value with the new estimate obtained using 'fminsearch':

    GaFrac1Resampled = interp1(t1, GaFrac1, tExpt, 'pchip');
    sse1 = norm(GaFracExpt - GaFrac1Resampled)^2;
    rSquare1 = 1 - sse1/sst
    
    rSquare1 =
    
        0.9199
  4. Plot the data and compare. If you left the previous figure open, because hold is on, the new plot appears in the existing figure to facilitate the comparison:

    plot(t1, GaFrac1, 'm-');
    legendText{end + 1} = sprintf('kGd Changed R^2=%4.2f', rSquare1);
    legend(legendText{:});

    The figure shows the best fit achieved by changing the parameter kGd.

      Note:   Leave this figure window open, so that you can use it later in this example.

Estimating Other Parameters in the G Protein Model

The example illustrating sensitivity analysis (Calculate Sensitivities) showed that Ga is sensitive to parameters kRs, kRD1, kGa, and kGd. Based on the results from the sensitivity analysis, this tutorial shows you how to estimate these parameters. The sensitivity data is presented in Extracting and Plotting Sensitivity Data.

    Note:   Although this example estimates four parameters to fit the data, there is no published experimental data that verifies these values, and this example is only for illustration.

  1. Create a variable containing the parameters to estimate:

    paramsToEst = [sbioselect(m1, 'Name', 'kRs');...
                   sbioselect(m1, 'Name', 'kRD1');...
                   sbioselect(m1, 'Name', 'kGa');...
                   sbioselect(m1, 'Name', 'kGd')];
    
  2. Estimate the parameters. Use the current values of parameters in the model as the starting values for optimization. Use the opt variable you created previously to specify plotting of each iteration of the parameter estimation to see how optimization is progressing:

    [estValues2, result2] = sbioparamestim(m1, tExpt, GaFracExpt,...
                            GaFrac, paramsToEst, {}, {'fminsearch',opt});
    

      Note:   Close this figure before proceeding with the example.

  3. Compare original parameter values and the estimated parameter values obtained with 'fminsearch':

    % Original parameter values
    paramsToEst
    
    SimBiology Parameter Array
    
    Index:    Name:    Value:    ValueUnits:
    1         kRs      4         
    2         kRD1     0.004     
    3         kGa      1e-005    
    4         kGd      0.11      
    
    % Estimated parameter values
    num2str(estValues2)
    
    ans =
    
          4.549
      0.0031018
    9.0068e-006
        0.12381
  4. Calculate the R2 value using the new estimates obtained with 'fminsearch':

    estvarObj2 = addvariant(m1, 'Optimized kRs, kRD1, kGa, and kGd');
    addcontent(estvarObj2, ...
        {{'parameter', 'kRs', 'Value', estValues2(1)}, ...
        {'parameter', 'kRD1', 'Value', estValues2(2)}, ...
        {'parameter', 'kGa', 'Value', estValues2(3)}, ...
        {'parameter', 'kGd', 'Value', estValues2(4)}});
    simDataObj2 = sbiosimulate(m1, estvarObj2);
    [t2, GaFrac2] = selectbyname(simDataObj2, 'GaFrac');
    GaFrac2Resampled = interp1(t2, GaFrac2, tExpt, 'pchip');
    sse2 = norm(GaFracExpt - GaFrac2Resampled)^2;
    rSquare2 = 1 - sse2/sst
    
    rSquare2 =
    
        0.9603
  5. Plot the data and compare. If you left the previous figure open, because hold is on, the new plot appears in the existing figure to facilitate the comparison:

    plot(t2, GaFrac2, 'g-');
    legendText{end + 1} = sprintf('4 Constants Changed R^2=%4.2f',...
                          rSquare2);
    legend(legendText{:});

Was this topic helpful?