How to fit a Hill model with the mosaic package?

Hello,
I have a set of data coming from a dissociation experiment
(protein/ligand). Since the data is required to calculate the constant
of dissociation (Kd) of this pair, I am looking for a way to fit a
Hill function to the data.
I have seen that the package `mosaic` (https://rpubs.com/dtkaplan/646) provides
such approach with `fitModel`.
However, I fail to find the proper starting points and get the error:
```
Error in nlsModel(formula, mf, start, wts, scaleOffset = scOff,
nDcentral = nDcntr) : singular gradient matrix at initial parameter
estimates
```
How can I run this package?
How can I estimate the starting values for the function, considering
that the data is decreasing rather than increasing (so it is probably
an inverse Hill function)?
Thank you

```
df = data.frame(Response =    c(890.72,    895.46,    895.63,
894.72,    895.49,    893.59,
                892.53,    895.06,    897.21,    889.27,    876.05,
857.96,    862.02,    858.36,
                890.94,    890.8,    887.22,    888.91,    890.83,
889.92,    891.76,    890.32,
                886.35,    878.11,    866.57,    859.04,    863.64,
880.16,    884.15,    879.57,
                878.89,    882.27,    881.59,    880.98,    881.45,
876.19,    868.32,    859.16,
                850.53,    853.21,    859.34,    859.73,    861.19),
                Dose =    c(0.0000000015,    0.000000003,
0.000000006,    0.000000012,
                       0.000000024,    0.000000048,    0.000000095,
0.00000018,
                       0.00000038,    0.00000078,    0.0000015,
0.000013,    0.000025,
                       0.00005,    0.0000000015,    0.000000003,
0.000000006,
                       0.000000012,    0.000000024,    0.000000048,
0.000000095,
                       0.00000018,    0.00000038,    0.00000078,
0.0000015,    0.000025,
                       0.00005,    0.0000000015,    0.000000003,
0.000000006,
                       0.000000012,    0.000000024,    0.000000048,
0.000000095,
                       0.00000018,    0.00000038,    0.00000078,
0.0000015,    0.000003,
                       0.000006,    0.000013,    0.000025,    0.00005)
)
plot(Response~log10(Dose), df)
nls(Response ~ ( (exp((a+b*Dose))) / (1+exp(-(a+b*Dose)) )), data=df)
library(mosaic)
fitModel(Response ~ A * Dose^2/(theta^2 + Dose^2), data = df,
         start = list(A = 100, theta = -4))
```