weighted non-linear fit help

I am trying to fit the following non-linear equation (SpO2) for two parameters S and HK1 where my input data is SpO2 (dependent range:[0.5,1]) and FIO2 (independent range:[0.5,1]). My S range is [0,1] (should normally be 0.1), my HK1 range is [0,+inf] (should normally be 1). I get a RootOf equation (caused by solving a cubic equation). If I substitute in values for FIO2, HK1 and S, I need to use the double command to have the answer solved. I will only need the real root, although 16 roots are returned, the correct one is #7 (verified elsewhere). I have tried to use fittype and nlinfit but they don't seem to want to take my equation as a "model". I am not terribly familiar with Matlab, so any help would be greatly appreciated.

Thanks!

syms FIO2 Hb L Pa PaCO2_1 Pb Pc1 Q Q_1 R RQ S Sa Sc1 T Va_1 a1 a2 a3 avD HK1
avD=5                       %known constant
Hb=15                       %known constant
a1=23400                    %known constant
a2=1                        %known constant
a3=150                      %known constant
Pb=760                      %known constant
RQ=0.8                      %known constant
L=863/100                   %known constant
Va_1 = HK1*Q_1              %A way to get a number for a ratio
eq1 = '134*Hb*Sc1*(1/100)+3*Pc1*(1/1000)-134*Hb*Sa*(1/100)-3*Pa*(1/1000) = S*avD/(1-S)'                               %Shunt equation
Sa = '1/(a1/(a2*Pa^3+a3*Pa)+1)'            %Severinghaus equation
Sc1 = '1/(a1/(a2*Pc1^3+a3*Pc1)+1)'
Pc1 = 'FIO2*(Pb-47)-PaCO2_1*(1-FIO2*(1-RQ))/RQ'
eq50a = 'Va_1/Q_1 = L*RQ*avD/PaCO2_1'
PaCO2_1=solve(eq50a,'PaCO2_1')
%Make the algebra easier to solve for Pa
Q='((Hb*Sc1*134)/(100))',R='((Pc1*3)/(1000))',T='(S/(1-S))*avD'
%Modified shunt equation
eq12='Q+R-(67/50)*Hb/(a1/(a2*Pa^3+a3*Pa)+1)-(3/1000)*Pa = T' 
Pa=solve(eq12,'Pa')                        %Rootof solution cause
Sap=subs(Sa,'Pa',Pa)
eqh1=subs(Sap,{'Q','R','T'},{Q,R,T})
%There is probably a more direct way to do this, substituting in all of my knowns and equations
eqh2=subs(eqh1,'Sc1',Sc1);
eqh2=subs(eqh2,'Pc1',Pc1)
eqh2=subs(eqh2,'PaCO2_1',PaCO2_1)
eqh2=subs(eqh2,'Va_1',Va_1)
SpO2=subs(eqh2,{'avD','Hb','a1','a2','a3','Pb','RQ','L'},{avD,Hb,a1,a2,a3,Pb,RQ,L})
%sample data set
FIO2data=[0.209 0.182 0.165 0.15 0.132 0.109 0.208 0.236 0.278 0.329]
SPO2data=[0.965 0.942 0.841 0.787 0.746 0.623 0.973 0.98 0.983 0.983]