Even in the context of fitting dose-response (or similar) curves, there are several types of “sigmoid” curves that might be fit. The “Alternative Hill Fit” shown in the second figure of the question is a special case of a more general log-logistic or Hill fit, with symmetry enforced.

The general form (in the terminology of the help page for LL.5 in the R drc package) is:

$$f(x)= c + \frac{d-c}{1+\exp(b(\log(x)-\log(e))^f},$$

where $c$ is the lower asymptote, $d$ is the upper asymptote, $e$ is related to a midpoint parameter, $b$ is related to the steepness of the curve between the asymptotes, and the exponent $f$ allows for asymmetry if its value is other than 1. The fit that’s shown can be matched very well by fixing $f=1$ first.

Manually recovering the mean values from the second figure of the question* gave:

meanValues
#  [1]  39.2  49.3  53.1  50.5  57.5  61.2  72.3  86.0 117.8 131.2
# [11] 139.8 138.7 132.0 131.4 130.4 128.2 128.8 125.6 125.9 125.1
meanCounts <- data.frame(mV=seq(100,2000,by=100),meanCount=meanValues)

Fitting the mean values to the above equation with $f=1$ came very close to the “Alternative Hill fit.”

library(drc)
fixedMeanHill <- drm(meanCount ~ mV, data=meanCounts,
                 fct=LL.5(fixed=c(NA,NA,NA,NA,1))) ## fixes f at 1
summary(fixedMeanHill)
# 
# Model fitted: Generalized log-logistic (ED50 as parameter) (4 parms)
# 
# Parameter estimates:
# 
#               Estimate Std. Error t-value   p-value  from "Alt Hill Fit"  
# b:(Intercept) -11.4432     2.6281 -4.3542 0.0004919  (not reported)
# c:(Intercept)  51.3728     2.8985 17.7238 6.111e-12  “Bottom = 50.86"
# d:(Intercept) 131.1763     1.9896 65.9295 < 2.2e-16  “Top = 131.50"
# e:(Intercept) 791.2281    21.0532 37.5823 < 2.2e-16  “EC50 = 794.38"
# 
# Residual standard error:
# 
#  6.305347 (16 degrees of freedom)

Allowing $f$ to vary, however, leads to a better fit.

fullMeanHill <- drm(meanCount ~ mV, data=meanCounts, fct=LL.5())
summary(fullMeanHill)
# several lines omitted
# Residual standard error:
# 
#  4.970377 (15 degrees of freedom)

ED(fullMeanHill,50) ## get ED50
# 
# Estimated effective doses
# 
#        Estimate Std. Error
# e:1:50  795.546     18.137

anova(fullMeanHill, fixedMeanHill)
##
# ANOVA table
# 
#           ModelDf    RSS Df F value p value
# 2nd model      16 636.12                   
# 1st model      15 370.57  1 10.7489  0.0051

This still doesn’t deal with three issues: the errors in the mean values, the apparent increase of standard errors with higher mean counts, and the apparent overshoot at mV levels around 1100 to 1200.

For the first two issues, counts of “spikes” in neural responses could be Poisson-distributed values (variance = mean). I generated 8 Poisson values around each of the mV values to get some variability for further modeling.

set.seed(20260721)
allCounts <- data.frame(mV = rep(seq(100, 2000, by = 100), each=8),
        counts = rpois(160,lambda = rep(meanCounts$meanCount,each = 8)))

The drc package can fit some generalized models. With a Poisson model for counts and more observations, the standard error of the ED50 estimate was lower:

fullPoissonHill <- drm(counts ~ mV, data = allCounts,
                       fct = LL.5(), type="Poisson")

ED(fullPoissonHill, 50)
# 
# Estimated effective doses
# 
#        Estimate Std. Error
# e:1:50  800.904     11.191

Sometimes there is an overshoot in a dose-response curve between the upper asymptote and the drop toward the lower asymptote. That’s called “hormesis” and can be modeled with a “Brain-Cousens” model. That adds a linear function of $x$ to the numerator of the general formula above while fixing the exponent $f=1$. For these synthetic data, at least, the fit is visually better.

fullBC <- drm(counts ~ mV, data = allCounts,
              fct = BC.5(),type = "Poisson")

fits of spike counts versus applied voltage

This isn’t classical hormesis, as this is a rising dose-response curve at values of mV below the overshoot (instead of a falling curve at values above the overshoot), and an ED50 thus isn’t available. The drop in counts at the highest levels of mV, if real, might be due to adaptation of the neural responses, particularly if the stimuli were always applied at increasing levels.

The drc package can fit several more forms of these curves, including logistic curves (with linear $x$ instead of log; not much different from log-logistic on these data) and classic Michaelis-Menten curves (which didn’t work well on these data).

Code for plot:

## get curves from these fits
p.FH <- plot(fixedMeanHill, log="")
names(p.FH) <- c("x","y")
p.PH <- plot(fullPoissonHill,log="")
names(p.PH) <- c("x","y")
p.BC <- plot(fullBC,log="")
names(p.BC) <- c("x","y")

library(ggplot2)
ggplot(data = allCounts, mapping = aes(x = mV, y = counts)) + 
stat_summary(fun = mean, geom = "point") + 
stat_summary(fun.data = mean_se, geom="errorbar", width=50) + 
geom_smooth(method ="gam", aes(col="gam smooth"), linewidth=1, linetype=2) + 
geom_line(data=p.FH, aes(x=x, y = y, col="Fixed Mean Hill"), linewidth=1.5) + 
geom_line(data=p.PH, aes(x=x, y = y, col="Poisson Hill"), linewidth=1.5) + 
geom_line(data=p.BC, aes(x=x, y = y, col="Brain-Cousens"), linewidth=2) +
scale_color_manual(name='Fit Type',
    breaks=c('gam smooth', 'Fixed Mean Hill', 'Poisson Hill', 'Brain-Cousens'),
    values=c('darkgray', 'yellow', 'blue', 'red'))

*The y-axis values for corresponding points seem to differ between the two figures in the question. AI-assisted data extraction from this plot didn’t work well for me.