a. Using only one of the MATLAB functions lsqcurvefit and fitnlm perform a nonlinear regression.

What are the optimal values for Emax and EC50 for each drug?

First, we need to write the data in the provided table as vectors in MATLAB syntax.

dose = [1 10 60 100 200 450 600 850 1000];

drug_A= [8 49 82 89 92 94 96 98 98];
drug_B= [8 31 58 75 77 78 80 79 84];

𝒄
Now we need to express the Hill model 𝑬(𝒄) = 𝑬𝒎𝒂𝒙 expression in MATLAB syntax.
𝑬𝑪𝟓𝟎 +𝒄

fitHillModel=@(b,x) b(1).*x./(b(2)+x);

where b(1) is Emax and b(2) is EC50.

Define the initial values to start with, they are just some random values.

init_values_A=[max(drug_A);median(dose)];

init_values_B=[max(drug_B);median(dose)];

mdl_A=lsqcurvefit(fitHillModel,init_values_A,dose,drug_A)

mdl_B=lsqcurvefit(fitHillModel,init_values_B,dose,drug_B)

By executing these lines of codes, we get the following optimal values for 𝑬𝒎𝒂𝒙 and 𝑬𝑪𝟓𝟎

Table 1: Values of Emax and EC50 for the non-linear regression

Drug 𝑬𝒎𝒂𝒙 𝑬𝑪𝟓𝟎

Drug A 97.6973 10.3186

Drug B 82.6374 17.5659

b. Linearize the Hill Equation using Lineweaver-Burk method. Using the linearized form of the
function and the data given above compute the EC50 value and the maximum effect (Emax)
for each drug. Are the values exactly the same as those you computed in (a)? Explain.

First let’s linearize mathematically. We take the reciprocal of E(C)

 1 𝐸𝐶50 + 𝑐 𝐸𝐶50 1 1
= = × +
𝐸(𝑐) 𝐸𝑚𝑎𝑥 × 𝑐 𝐸𝑚𝑎𝑥 𝑐 𝐸𝑚𝑎𝑥

Where B0 is 1/Emax and B1 is EC50/Emax.

Therefore, Emax = 1/B0 and EC50= Emax * B1

Check the code in the appendix for the full implementation. Results are presented in the

following table.

Table 2: Values of Emax and EC50 for the linear regression

linear regression E max EC50

Drug_A 1/ 0.010058 = 99.4233446 99.4233446 * 0.11484 =11.41777689

Drug_B 1/ 0.013692=73.03534911 73.03534911 * 0.11204 =8.182880514

By Comparing the results in (a) and (b), it is obvious that there is a difference between the linear

and the non-linear approach since linear regression is just an approximation while the hill model

is a representation of the exact function. Linearization thus will make us lose some of the

information. Thus, the data in (a) are more accurate and reliable.
 c. Plot separately the fitted curve and the fitted line that you computed in (a) and (b),

together with the data used. What do you observe?

Figure 1: Drug A Hill model and the linear regression fit.

Figure 2: Drug B Hill model and the linear regression fit.

In figure 1, initially we do observe that as we increase the dose concentration, the effect of the
drug also increases. However, at some point the curve reaches a stationary phase. These findings do
make sense since drugs have a certain maximum effect. In drug A for example, we see that once the
dose concentration exceeds 400, the effect of the drug stays the same. This aligns with what we have
found in part a since we found Emax value to be ~97. In other words, despite the dose provided, the
maximum effect that drug A can have is 97 and that is obvious in the plot.
 In addition, by comparing the hill models in figures 1 and 2, we see that drug has a higher
maximum effect and seems to be better than drug B.

In the linear regression fit of the same drug A, we observe that there’s a linear relationship
between the dose concentration and the effect of the drug. As we increase the dose from 0 to 1000, the
effect of the drug increases from 0 to approximately 97. However, we know that this is not the situation,
and the non-linear model is more reliable as we proved in part b.

The same analysis logic applies for drug B.

d. Compute 𝑹𝟐 for both (a) and (b). Which fit seems to be better?

Already, we got the 𝑹𝟐 values for the linear models, and the implementation of the code in order to
get the 𝑹𝟐 for the non-linear models is in the appendix under part d. The 𝑹𝟐 values for the models are
presented in the following table.
Table 3: Summary of 𝑹𝟐 values

Models Linear Regression Non-linear Linear regression Non-linear

(A) Regression (A) (B) Regression (B)
𝑹𝟐 1 0.998 0.993 0.983

As a rule of thumb, the higher the 𝑹𝟐 value, the better model we have. We do observe that the
𝟐
𝑹 values are being higher in the linear regression model for both drugs (A being the highest). However,
when it comes to interpreting these results, it is better to look at the 𝑹𝟐 values in the non-linear models
since these models are more reliable. Having said that, it is concluded that the model for drug A fits
better.

e. Considering the definitions below, compare the efficacy and potency of the drugs. Which one
would you prefer and why?

Efficacy is related to the maximum effect of the drug regardless of the dose provided. Thus, the higher
the efficacy (Emax) of a drug, the better the drug.

Potency is the required concentration of a drug to produce a given effect, the potency of a drug is
related to EC50.

Therefore, the best drug would be the one with high effect (𝐸𝑚𝑎𝑥 ) and low 𝐸𝐶50

By looking at the values of 𝐸𝑚𝑎𝑥 and 𝐸𝐶50 summarized in table 1 and 2. We conclude that in both
models, drug A seems to be the most efficient and potent drug.
1 APPENDIX

1.1 CODE

% PART a
%*****************
%entering the data provided in the table.

dose = [1,10,60,100,200,450,600,850,1000];
drug_A= [8 49 82 89 92 94 96 98 98];
drug_B= [8 31 58 75 77 78 80 79 84];
%expressing the hill model in MATLAB syntax
fitHillModel=@(b,x) b(1).*x./(b(2)+x);
%Defining some initial values to start with, just some random values \

init_values_A=[max(drug_A);median(dose)];
init_values_B=[max(drug_B);median(dose)];

mdl_A=lsqcurvefit(fitHillModel,init_values_A,dose,drug_A)
mdl_B=lsqcurvefit(fitHillModel,init_values_B,dose,drug_B)
mdl_A =

97.6973
10.3186

mdl_B =

82.6374
17.5659

%PART b
%*******************
%Linearize the equation by taking the reciprocal of the dose

dose_linear= 1./dose';
drug_A_linear=1./drug_A';
drug_B_linear=1./drug_B';
%put them in table
table=table(dose_linear,drug_A_linear,drug_B_linear);
%The A linear model
mdl_A_linear=fitlm(table,'drug_A_linear~dose_linear')

mdl_A_linear =

Linear regression model:

drug_A_linear ~ 1 + dose_linear

Estimated Coefficients:
Estimate SE tStat pValue
 ________ __________ ______ __________

(Intercept) 0.010058 0.0001687 59.62 9.8035e-11

dose_linear 0.11484 0.00050348 228.08 8.2217e-15

Number of observations: 9, Error degrees of freedom: 7

Root Mean Squared Error: 0.000469
R-squared: 1, Adjusted R-Squared: 1
F-statistic vs. constant model: 5.2e+04, p-value = 8.22e-15
mdl_B_linear=fitlm(table,'drug_B_linear~dose_linear')

mdl_B_linear =

Linear regression model:

drug_B_linear ~ 1 + dose_linear

Estimated Coefficients:
Estimate SE tStat pValue
________ _________ ______ __________

(Intercept) 0.013692 0.0011339 12.075 6.0986e-06

dose_linear 0.11204 0.0033842 33.107 5.9389e-09

Number of observations: 9, Error degrees of freedom: 7

Root Mean Squared Error: 0.00315
R-squared: 0.994, Adjusted R-Squared: 0.993
F-statistic vs. constant model: 1.1e+03, p-value = 5.94e-09

%PART c
%******************
%a vector with range [minumum dose, maximum dose]
AgVec=linspace(min(dose),max(dose));
%The pridicted respone
predResponse=fitHillModel(mdl_A,AgVec);
%create a figure
fig1=figure(1);
%Plot the response of drug A Vs dose (stars, size of 10)
plot(dose,drug_A,'bp','MarkerSize', 10)
%hold on and keep the changes on the same graph
hold on
%Plot the curve
plot(AgVec,predResponse,'-r')

%linearized form of model A. Create a vector with values of Coefficients

coeffs = mdl_A_linear.Coefficients.Estimate;
%X vector ranging from min dose to max dose
x1fit = min(dose):1:max(dose);
y1fit = coeffs(1)+coeffs(2)*x1fit;
%Plot the model
plot(x1fit,y1fit,'b-')
%add lables on the graph.
legend('Data','Hill Equation Fit','Linear fit','Location','SE')
xlabel('Dose')
ylabel('Response of Drug A')
% we repeat the same for Drug B only changing the names of the varaibles.

%PART d
%*********************
%agin some initial values to start with
initial_A=[max(drug_A);median(dose)];
%use the fitlm function to compute the R^2 value for the non-linear model
nonlinearA= fitnlm(dose,drug_A,fitHillModel,initial_A)

nonlinearA =

Nonlinear regression model:

y ~ b1*x/(b2 + x)

Estimated Coefficients:
Estimate SE tStat pValue
________ _______ ______ __________

b1 97.697 0.55531 175.93 5.0586e-14

b2 10.319 0.4882 21.136 1.3353e-07

Number of observations: 9, Error degrees of freedom: 7

Root Mean Squared Error: 1.21
R-Squared: 0.999, Adjusted R-Squared 0.998
F-statistic vs. zero model: 2.15e+04, p-value = 5.53e-14
initial_B=[max(drug_B);median(dose)];
nonlinearB= fitnlm(dose,drug_B,fitHillModel,initial_B)

nonlinearB =

Nonlinear regression model:

y ~ b1*x/(b2 + x)

Estimated Coefficients:
Estimate SE tStat pValue
________ ______ ______ __________

b1 82.638 1.7795 46.44 5.6128e-10

b2 17.566 2.8018 6.2697 0.00041617

Number of observations: 9, Error degrees of freedom: 7

Root Mean Squared Error: 3.51
R-Squared: 0.985, Adjusted R-Squared 0.983
F-statistic vs. zero model: 1.69e+03, p-value = 3.98e-10