See lsqcurvefit(): https://www.mathworks.com/help/releases/R2020a/optim/ug/lsqcurvefit.html. You can specify lb and ub parameters to specify the lower and upper bounds.
For example
x = [6.1,6.1,6.1,6.1,6.1,6.1,6.1,6.4,6.4,6.4,6.4,6.4,6.4,6.4,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7,6.7];
y = [10,20,30,40,50,60,70,10,20,30,40,50,60,70,10,20,30,40,50,60,70,10,20,30,40,50,60,70];
z = [0.142971,0.07939,0.051199,0.040022,0.033678,0.041687,0.025369,0.376911,0.255986,0.222872,0.178438,0.140919,0.11479,0.09615,0.310793,0.190409,0.183967,0.132389,0.086559,0.083365,0.077659,0.932289,0.466257,0.282995,0.156609,0.145649,0.126456,0.126298];
% p is defined as [a b c d e]
f = @(p, xdata) p(1)*(xdata(:,1)+p(2))+p(3)*log10(xdata(:,2)+p(4))+p(5);
xdata = [x(:) y(:)];
ydata = z(:);
p_lb = [0 -10 -5 -5 0]; % lower bounds
p_ub = [1 0 0 0 1]; % upper bounds
p_sol = lsqcurvefit(f, zeros(1,5), xdata, ydata, p_lb, p_ub);
z_est = f(p_sol, xdata);
plot(z);
hold on
plot(z_est);
