Fitting data to a Boltzmann sigmoid function, and estimating parameters with a non-linear least squares fitting in R

I have a dataset of 10 subjects, whose muscle response amplitude was measured across 11 different stimulus intensities (states) applied to their muscles. The output curve of the data looks like a sigmoidal function, with different stimulus intensities on the x-axis and resultant muscle amplitude outputs on the y-axis.

I would like to fit a Boltzmann sigmoid function to fit each subject's data. From the literature in my field, I would need to estimate my parameters of interest (i.e., slope, S50, and plateau from the function below) using a non-linear least squares fitting using R (typical once used is a Marquardt-Levenberg algorithm for least squares convergence).

From what I could find, I believe the Boltzmann function should look like:

y(x) = plateau / [1 + e^((S50-x)/slope)]

Plateau = maximum amplitude
S50 = stimulus intensity x required to evoke a response equal to half the maximum amplitude

Example data for one subject:

y = c(0.04, 0.053, 0.096, 0.219, 0.701, 1.333, 3.032, 3.556,4.33, 4.27, 4.167)
x = c(1,2,3,4,5,6,7,8,9,10,11)

Giving:

This post has someone outline how they do this with the "nls" package: https://stackoverflow.com/questions/60800267/working-with-nls-function-produces-a-singular-gradient/60800647#60800647

So I would have something like:

fit <- nls(y ~ plateau / (1 + exp((S50 - x)/slope), start = list( ??? ))

However, I'm still a little lost. How do I know what parameters to feed into the equation (which I think would go in the [???] section)? Do I need to specify plateau, S50, and slope? Isn't the whole point that I don't know these going in?

Is there a way to feed my data into some function to then get an output of the best non-linear equation that fits my data and gives me those S50, plateau, and slope parameters? Am I misunderstanding the process?

There is apparently a "gslnls" package available in R that uses the Marquardt-Levenberg algorithm to fit the data. I guess I am a little confused as to how to apply it, or if there is a simpler way to do this.

Any guidance would be greatly appreciated.

What goes inside the [???] are the initial guesses for the starting values for your three parameters. R and nls() does not know how your function behaves nor where are feasible solutions exist thus thus it needs some help on where to start looking. I recommend adding the trace=TRUE option to function to see where is it looking for the solution. Sometimes this takes some trial and error to get right!

Singular gradient error usually means one is entering an undefined region of the solution space, i.e. a value is going to infinity, square root of a negative number, the function is over defined with the number of adjustable parameters, etc.

Starting with the slope <1 seems to work.

y = c(0.04, 0.053, 0.096, 0.219, 0.701, 1.333, 3.032, 3.556,4.33, 4.27, 4.167)
x = c(1,2,3,4,5,6,7,8,9,10,11)
fit <- nls(y ~ plateau / (1 + exp((S50 - x)/slope)), start = list( plateau=1, S50=1, slope=.5), trace = TRUE)