Legend with ggplot and geom_sf

Question

I have built a map of Toulouse using ggplot and geom_sf for having only the elements I want, eg. main roads, river, railways, ... Then I need to plot the n bike stations, each one colored with the n values of a vector F0. These colors constitute a gradient from white to darkred. But I don't arrive to plot the gradient of colors in a legend. If someone can help!

My code is the following:

library(osmdata) # wrapper for Overpass API from Open Street Maps
library(sf) # library for manipulating Simple features objects
library(dplyr) 
library(ggplot2)
toulouse <- opq("Toulouse")
head(toulouse)

#routes
toulouse_street1 <-  toulouse %>%
  add_osm_feature(key = "highway", 
                  value = c("motorway","motorway_link")) %>%
  osmdata_sf()
#### water
toulouse_water1<- toulouse %>%
  add_osm_feature(key = "natural", 
                  value = c("water")) %>%
  osmdata_sf()
toulouse_water1$osm_polygons <- toulouse_water1$osm_polygons %>% 
  filter(natural == "water")

all_toulouse_water <- st_union(toulouse_water1$osm_polygons,toulouse_water1$osm_multipolygons)

#### train et métro
toulouse_railway1 <- toulouse %>%
  add_osm_feature(key = "public_transport", value=c("station")) %>%
  osmdata_sf()
toulouse_railway2 <- toulouse %>%
  add_osm_feature(key = "railway",value=c("rail")) %>%
  osmdata_sf()

#####
#  loadings
######

library(sf)
library(autoimage)
shp <- st_read(dsn = '/Users',
               layer = 'factor_loadings_with_geometry_T2')
coord<-shp$geometry
coord_new<-st_transform(coord,crs=4326)
F0<-shp$F5 #F0 is between 0 and 1
mycol<-colorRampPalette(c("white","bisque","indianred1","red3","red4","grey20"))
coul=mycol(11)
F0_ok=numeric()
F0_ok=coul[F0+1]

###########

### print the plot
street_plot <- 
  ggplot() +
  geom_sf(data = toulouse_railway2$osm_lines,
          inherit.aes = FALSE,
          size = 0.2, 
          colour = "navajowhite3")+
  geom_sf(data = toulouse_railway2$osm_lines,
          inherit.aes = FALSE,
          size = 0.2, 
          colour = "grey35",
          linetype="dotted")+
  geom_sf(data = toulouse_street1$osm_lines,
          inherit.aes = FALSE,
          color = "grey40",
          size = 1) +
  geom_sf(data = toulouse_water1$osm_multipolygons,
          inherit.aes = FALSE,
          fill="skyblue1",
          alpha = .8,
          lwd=0) +
  geom_sf(data = coord_new,
          color = F0_ok,#rgb(1-F0, F0*F0, F0*F0, maxColorValue = 1), 
          fill = F0_ok,#rgb(1-F0, F0*F0, F0*F0, maxColorValue = 1), 
          shape = 21,
          size = 1.4,
          inherit.aes = FALSE
  )+
  annotation_north_arrow(which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"),location="bl",pad_x = unit(0.25, "cm"), pad_y = unit(0.5, "cm")) +
  annotation_scale(height = unit(0.15, "cm"))+
  #coord Toulouse
  coord_sf(xlim = c(1.39, 1.495758), ylim = c(43.55, 43.65)) +
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.line = element_blank(),
        panel.grid=element_blank(),
        panel.background = element_rect(fill = "floralwhite"),
        legend.position = "right")
# Print the plot
street_plot