diff --git a/NAMESPACE b/NAMESPACE
index 34b8a6e2..6d27aabb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -146,6 +146,7 @@ export(setXGrid)
 export(setYAxis)
 export(setYGrid)
 export(tlfStatFunctions)
+export(updateTimeProfileLegend)
 export(useDarkTheme)
 export(useExcelTheme)
 export(useHighChartTheme)
diff --git a/R/aaa-utilities.R b/R/aaa-utilities.R
index b05fce92..6e15cdff 100644
--- a/R/aaa-utilities.R
+++ b/R/aaa-utilities.R
@@ -87,10 +87,10 @@ parseUpdateAestheticProperty <- function(aestheticProperty, plotConfigurationPro
     parse(text = paste0(aestheticProperty, "Length <- length(unique(mapData[, ", aestheticProperty, "Variable]))")),
     # Update the property using ggplot `scale` functions
     parse(text = paste0(
-      "plotObject <- plotObject + ggplot2::scale_", aestheticProperty, "_manual(",
+      "suppressMessages(plotObject <- plotObject + ggplot2::scale_", aestheticProperty, "_manual(",
       "values=getAestheticValues(n=", aestheticProperty, "Length,",
       "selectionKey=plotConfiguration$", plotConfigurationProperty, "$", aestheticProperty,
-      ',aesthetic = "', aestheticProperty, '"))'
+      ',aesthetic = "', aestheticProperty, '")))'
     )),
     # remove the legend of aesthetic if default unmapped aesthetic
     parse(text = paste0("if(isIncluded(", aestheticProperty, 'Variable, "legendLabels")){plotObject <- plotObject + ggplot2::guides(', aestheticProperty, " = 'none')}"))
diff --git a/R/datamapping-range.R b/R/datamapping-range.R
index 95661108..d7d0e828 100644
--- a/R/datamapping-range.R
+++ b/R/datamapping-range.R
@@ -21,6 +21,7 @@ RangeDataMapping <- R6::R6Class(
     #' @param linetype R6 class `Grouping` object or its input
     #' @param shape R6 class `Grouping` object or its input
     #' @param size R6 class `Grouping` object or its input
+    #' @param group R6 class `Grouping` object or its input
     #' @param data data.frame to map used by `smartMapping`
     #' @return A new `RangeDataMapping` object
     initialize = function(x = NULL,
@@ -32,6 +33,7 @@ RangeDataMapping <- R6::R6Class(
                           linetype = NULL,
                           shape = NULL,
                           size = NULL,
+                          group = NULL,
                           data = NULL) {
 
       # smartMapping is available in utilities-mapping.R
@@ -39,7 +41,7 @@ RangeDataMapping <- R6::R6Class(
       super$initialize(x %||% smartMap$x,
         y = NULL,
         groupMapping = groupMapping, color = color, fill = fill,
-        linetype = linetype, shape = shape, size = size
+        linetype = linetype, shape = shape, size = size, group = group
       )
 
       self$ymin <- ymin %||% smartMap$ymin
diff --git a/R/datamapping-xygroup.R b/R/datamapping-xygroup.R
index 5a8359cc..ec01a449 100644
--- a/R/datamapping-xygroup.R
+++ b/R/datamapping-xygroup.R
@@ -18,6 +18,7 @@ XYGDataMapping <- R6::R6Class(
     #' @param linetype R6 class `Grouping` object or its input
     #' @param shape R6 class `Grouping` object or its input
     #' @param size R6 class `Grouping` object or its input
+    #' @param group R6 class `Grouping` object or its input
     #' @param data data.frame to map used by `smartMapping`
     #' @return A new `XYGDataMapping` object
     initialize = function(x = NULL,
@@ -28,6 +29,7 @@ XYGDataMapping <- R6::R6Class(
                           linetype = NULL,
                           shape = NULL,
                           size = NULL,
+                          group = NULL,
                           data = NULL) {
 
       # smartMapping is available in utilities-mapping.R
@@ -39,15 +41,16 @@ XYGDataMapping <- R6::R6Class(
         "fill" = fill,
         "linetype" = linetype,
         "shape" = shape,
-        "size" = size
+        "size" = size,
+        "group" = group
       ))
       # To simplify the process workflow, groupMapping inputs color, fill... can be used directly instead of groupMapping
       self$groupMapping <- groupMapping %||% GroupMapping$new(
-        color %||% smartMap$color,
-        fill %||% smartMap$fill,
-        linetype %||% smartMap$linetype,
-        shape %||% smartMap$shape,
-        size %||% smartMap$size
+        color %||% group %||% smartMap$color,
+        fill %||% group %||% smartMap$fill,
+        linetype %||% group %||% smartMap$linetype,
+        shape %||% group %||% smartMap$shape,
+        size %||% group %||% smartMap$size
       )
     },
 
diff --git a/R/observed-data-mapping.R b/R/observed-data-mapping.R
index b5cf7345..049b6758 100644
--- a/R/observed-data-mapping.R
+++ b/R/observed-data-mapping.R
@@ -6,8 +6,6 @@ ObservedDataMapping <- R6::R6Class(
   "ObservedDataMapping",
   inherit = XYGDataMapping,
   public = list(
-    #' @field lloq mapping lower limit of quantitation variable
-    lloq = NULL,
     #' @field error mapping error bars around scatter points
     error = NULL,
     #' @field mdv mapping missing dependent variable
@@ -21,8 +19,9 @@ ObservedDataMapping <- R6::R6Class(
     #' @param x Name of x variable to map
     #' @param y Name of y variable to map
     #' @param group R6 class `Grouping` object or its input
+    #' @param color R6 class `Grouping` object or its input
+    #' @param shape R6 class `Grouping` object or its input
     #' @param data data.frame to map used by `smartMapping`
-    #' @param lloq mapping lower limit of quantitation variable
     #' @param uncertainty mapping error bars around scatter points.
     #' Deprecated parameter replaced by `error`.
     #' @param error mapping error bars around scatter points
@@ -36,19 +35,18 @@ ObservedDataMapping <- R6::R6Class(
                           error = NULL,
                           ymin = NULL,
                           ymax = NULL,
-                          lloq = NULL,
                           mdv = NULL,
+                          color = NULL,
+                          shape = NULL,
                           group = NULL,
                           data = NULL){
-      validateIsString(lloq, nullAllowed = TRUE)
       validateIsString(uncertainty, nullAllowed = TRUE)
       validateIsString(error, nullAllowed = TRUE)
       validateIsString(ymin, nullAllowed = TRUE)
       validateIsString(ymax, nullAllowed = TRUE)
       validateIsString(mdv, nullAllowed = TRUE)
-      super$initialize(x = x, y = y, color = group, data = data)
+      super$initialize(x = x, y = y, color = color, shape = shape, group = group, data = data)
       
-      self$lloq <- lloq
       # If defined, ymin and ymax are used as is
       # If not, error/uncertainty are used and 
       # creates ymin and ymax as y +/- error
@@ -66,19 +64,14 @@ ObservedDataMapping <- R6::R6Class(
     checkMapData = function(data, metaData = NULL) {
       validateIsOfType(data, "data.frame")
       validateIsIncluded(self$error, names(data), nullAllowed = TRUE)
-      validateIsIncluded(self$lloq, names(data), nullAllowed = TRUE)
       validateIsIncluded(self$mdv, names(data), nullAllowed = TRUE)
       
       # Using super method, fetches x, y and groups
       mapData <- super$checkMapData(data, metaData)
-      # Add lloq data
-      if (!isOfLength(self$lloq, 0)) {
-        mapData[, self$lloq] <- data[, self$lloq]
-        mapData$lloq <- data[, self$lloq]
-      }
+      
       # ymin and ymax for error bars
       # This section may change depending on how we want to include options
-      if (!isOfLength(self$error, 0)) {
+      if (!isEmpty(self$error)) {
         mapData[, self$error] <- data[, self$error]
         # Symetric error bars
         mapData[, self$ymax] <- data[, self$y] + data[, self$error]
@@ -92,7 +85,7 @@ ObservedDataMapping <- R6::R6Class(
         mapData[, self$ymin] <- data[,self$ymin]
       }
       # MDV is a Nonmem notation in which values with MDV==1 are removed
-      if (!isOfLength(self$mdv, 0)) {
+      if (!isEmpty(self$mdv)) {
         mapData[, self$mdv] <- as.logical(data[,self$mdv])
         mapData <- mapData[!mapData[, self$mdv], ]
       }
diff --git a/R/plot-timeprofile.R b/R/plot-timeprofile.R
index 4550af91..181d763b 100644
--- a/R/plot-timeprofile.R
+++ b/R/plot-timeprofile.R
@@ -3,37 +3,36 @@
 #' Producing Time Profile plots
 #'
 #' @inheritParams addScatter
-#' @param dataMapping 
+#' @param dataMapping
 #' A `TimeProfileDataMapping` object mapping `x`, `y`, `ymin`, `ymax` and aesthetic groups to their variable names of `data`.
 #' @param observedData A data.frame to use for plot.
 #' Unlike `data`, meant for simulated data, plotted as lines and ribbons;
 #' `observedData` is plotted as scatter points and errorbars.
-#' @param observedDataMapping 
-#' An `ObservedDataMapping` object mapping `x`, `y`, `ymin`, `ymax`, `lloq` and aesthetic groups to their variable names of `observedData`.
-#' @param plotConfiguration 
+#' @param observedDataMapping
+#' An `ObservedDataMapping` object mapping `x`, `y`, `ymin`, `ymax` and aesthetic groups to their variable names of `observedData`.
+#' @param plotConfiguration
 #' An optional `TimeProfilePlotConfiguration` object defining labels, grid, background and watermark.
 #' @return A `ggplot` object
 #'
 #' @export
 #' @family molecule plots
-#' @examples 
+#' @examples
 #' # Produce a Time profile plot with observed and simulated data
 #' obsData <- data.frame(x = c(1, 2, 1, 2, 3), y = c(5, 0.2, 2, 3, 4))
 #' simTime <- seq(1, 10, 0.1)
 #' simData <- data.frame(
-#' x = simTime, 
-#' y = 10*exp(-simTime),
-#' ymin = 8*exp(-simTime),
-#' ymax = 12*exp(-simTime)
+#'   x = simTime,
+#'   y = 10 * exp(-simTime),
+#'   ymin = 8 * exp(-simTime),
+#'   ymax = 12 * exp(-simTime)
 #' )
-#' 
+#'
 #' plotTimeProfile(
-#' data = simData, 
-#' observedData = obsData,
-#' dataMapping = TimeProfileDataMapping$new(x = "x", y = "y", ymin = "ymin", ymax = "ymax"),
-#' observedDataMapping = ObservedDataMapping$new(x = "x", y = "y")
+#'   data = simData,
+#'   observedData = obsData,
+#'   dataMapping = TimeProfileDataMapping$new(x = "x", y = "y", ymin = "ymin", ymax = "ymax"),
+#'   observedDataMapping = ObservedDataMapping$new(x = "x", y = "y")
 #' )
-#' 
 plotTimeProfile <- function(data = NULL,
                             metaData = NULL,
                             dataMapping = NULL,
@@ -43,15 +42,15 @@ plotTimeProfile <- function(data = NULL,
                             plotObject = NULL) {
   validateIsOfType(data, "data.frame", nullAllowed = TRUE)
   validateIsOfType(observedData, "data.frame", nullAllowed = TRUE)
-  if (all(isOfLength(data, 0), isOfLength(observedData, 0))) {
-    warning("data and observed data are of length 0. Time profile layer was not added.")
+  if (all(isEmpty(data), isEmpty(observedData))) {
+    warning("'data' and 'observedData' are of length 0. Time profile layer was not added.")
     return(plotObject)
   }
 
-  if (!isOfLength(data, 0)) {
+  if (!isEmpty(data)) {
     dataMapping <- dataMapping %||% TimeProfileDataMapping$new(data = data)
   }
-  if (!isOfLength(observedData, 0)) {
+  if (!isEmpty(observedData)) {
     observedDataMapping <- observedDataMapping %||% ObservedDataMapping$new(data = data)
   }
 
@@ -65,60 +64,210 @@ plotTimeProfile <- function(data = NULL,
   plotObject <- plotObject %||% initializePlot(plotConfiguration)
 
   # Get transformed data from mapping and convert labels into characters usable by aes_string
-  if (!isOfLength(data, 0)) {
+  if (!isEmpty(data)) {
     mapData <- dataMapping$checkMapData(data)
-    if (!any(isOfLength(dataMapping$ymin, 0), isOfLength(dataMapping$ymax, 0))) {
-      plotObject <- addRibbon(
-        data = mapData,
-        dataMapping = dataMapping,
-        plotConfiguration = plotConfiguration,
-        plotObject = plotObject
-      )
+    mapLabels <- getAesStringMapping(dataMapping)
+    if (!any(isEmpty(dataMapping$ymin), isEmpty(dataMapping$ymax))) {
+      plotObject <- plotObject +
+        ggplot2::geom_ribbon(
+          data = mapData,
+          mapping = ggplot2::aes_string(
+            x = mapLabels$x,
+            ymin = mapLabels$ymin,
+            ymax = mapLabels$ymax,
+            fill = mapLabels$fill
+          ),
+          alpha = getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$alpha, position = 0, aesthetic = "alpha"),
+          show.legend = TRUE
+        )
     }
-    if (!isOfLength(dataMapping$y, 0)) {
-      plotObject <- addLine(
-        data = mapData,
-        dataMapping = dataMapping,
-        plotConfiguration = plotConfiguration,
-        plotObject = plotObject
-      )
+    if (!isEmpty(dataMapping$y)) {
+      plotObject <- plotObject +
+        ggplot2::geom_path(
+          data = mapData,
+          mapping = ggplot2::aes_string(
+            x = mapLabels$x,
+            y = mapLabels$y,
+            color = mapLabels$color,
+            linetype = mapLabels$linetype
+          ),
+          size = getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = 0, aesthetic = "size"),
+          show.legend = TRUE
+        )
     }
+    eval(parseUpdateAestheticProperty(AestheticProperties$fill, "ribbons"))
+    eval(parseUpdateAestheticProperty(AestheticProperties$linetype, "lines"))
+  }
+
+  # If no observed data, also update colors and return plotObect
+  if (isEmpty(observedData)) {
+    eval(parseUpdateAestheticProperty(AestheticProperties$color, "lines"))
+    plotObject <- setLegendPosition(plotObject)
+    plotObject <- setLegendFont(plotObject)
+    eval(parseUpdateAxes())
+    return(plotObject)
   }
 
-  if (!isOfLength(observedData, 0)) {
-    mapObservedData <- observedDataMapping$checkMapData(observedData)
-    if (!all(
-      isOfLength(observedDataMapping$error, 0), 
-      isOfLength(observedDataMapping$ymin, 0),
-      isOfLength(observedDataMapping$ymax, 0)
-    )){
-      plotObject <- addErrorbar(
+  mapObservedData <- observedDataMapping$checkMapData(observedData)
+  observedMapLabels <- getAesStringMapping(observedDataMapping)
+
+  if (!any(isEmpty(observedDataMapping$ymin), isEmpty(observedDataMapping$ymax))) {
+    # Split errorbars for negative data and log scaling
+    plotObject <- plotObject +
+      ggplot2::geom_linerange(
         data = mapObservedData,
-        dataMapping = observedDataMapping,
-        plotConfiguration = plotConfiguration,
-        plotObject = plotObject
-      )
-    }
-    plotObject <- addScatter(
-      data = mapObservedData,
-      dataMapping = observedDataMapping,
-      plotConfiguration = plotConfiguration,
-      plotObject = plotObject
-    )
-    # LLOQ
-    if (!isOfLength(observedDataMapping$lloq, 0)) {
-      mapObservedData$legendLabels <- paste(mapObservedData$legendLabels, " LLOQ")
-      plotObject <- addLine(
+        mapping = ggplot2::aes_string(
+          x = observedMapLabels$x,
+          ymin = observedMapLabels$ymin,
+          ymax = observedMapLabels$y,
+          color = observedMapLabels$color
+        ),
+        size = getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$size, position = 0, aesthetic = "size"),
+        linetype = getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$linetype, position = 0, aesthetic = "linetype"),
+        show.legend = FALSE
+      ) +
+      ggplot2::geom_linerange(
         data = mapObservedData,
-        dataMapping = XYGDataMapping$new(
-          x = observedDataMapping$x,
-          y = "lloq",
-          color = "legendLabels"
+        mapping = ggplot2::aes_string(
+          x = observedMapLabels$x,
+          ymin = observedMapLabels$y,
+          ymax = observedMapLabels$ymax,
+          color = observedMapLabels$color
         ),
-        plotConfiguration = plotConfiguration,
-        plotObject = plotObject
+        size = getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$size, position = 0, aesthetic = "size"),
+        linetype = getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$linetype, position = 0, aesthetic = "linetype"),
+        show.legend = FALSE
       )
+  }
+  plotObject <- plotObject +
+    ggplot2::geom_point(
+      data = mapObservedData,
+      mapping = ggplot2::aes_string(
+        x = observedMapLabels$x,
+        y = observedMapLabels$y,
+        color = observedMapLabels$color,
+        shape = observedMapLabels$shape
+      ),
+      size = getAestheticValues(n = 1, selectionKey = plotConfiguration$points$size, position = 0, aesthetic = "size"),
+      show.legend = TRUE
+    )
+
+  # Update shapes
+  # Code chunk below is equivalent to commented expression with a change of variable names
+  # parseUpdateAestheticProperty(AestheticProperties$shape, "points")
+  shapeVariable <- gsub("`", "", observedMapLabels$shape)
+  shapeLength <- length(levels(mapObservedData[, shapeVariable]))
+  shapeValues <- getAestheticValues(
+    n = shapeLength,
+    selectionKey = plotConfiguration$points$shape,
+    aesthetic = "shape"
+  )
+  suppressMessages(plotObject <- plotObject + ggplot2::scale_shape_manual(values = shapeValues))
+  if (isIncluded(shapeVariable, "legendLabels")) {
+    plotObject <- plotObject + ggplot2::guides(shape = "none")
+  }
+
+  # Update colors,
+  # Since colors can be available in both simulated and observed, the commented expressions can't apply
+  # parseUpdateAestheticProperty(AestheticProperties$color, "lines")
+  # parseUpdateAestheticProperty(AestheticProperties$color, "points")
+
+  # No simulated data -> update only observedData
+  if (isEmpty(data)) {
+    colorVariable <- gsub("`", "", observedMapLabels$color)
+    colorLength <- length(levels(mapObservedData[, colorVariable]))
+
+    suppressMessages(
+      plotObject <- plotObject +
+        ggplot2::scale_color_manual(
+          values = getAestheticValues(
+            n = colorLength,
+            selectionKey = plotConfiguration$points$color,
+            aesthetic = "color"
+          )
+        )
+    )
+    if (isIncluded(colorVariable, "legendLabels")) {
+      plotObject <- plotObject + ggplot2::guides(color = "none")
     }
+
+    plotObject <- setLegendPosition(plotObject)
+    plotObject <- setLegendFont(plotObject)
+    eval(parseUpdateAxes())
+    return(plotObject)
+  }
+
+  # Simulated and Observed data -> need to merge the legends
+  colorVariable <- gsub("`", "", mapLabels$color)
+  colorLength <- length(levels(mapData[, colorVariable]))
+  colorObservedVariable <- gsub("`", "", observedMapLabels$color)
+  # The final color vector needs a length of totalLength to prevent scale_color_manual to crash
+  colorBreaks <- c(
+    levels(mapData[, colorVariable]), 
+    setdiff(levels(mapObservedData[, colorObservedVariable]), levels(mapData[, colorVariable]))
+  )
+  totalLength <- length(colorBreaks)
+
+  # colorValues is colors for simulated and then colors for observed
+  colorValues <- c(
+    getAestheticValues(
+      n = colorLength,
+      selectionKey = plotConfiguration$lines$color,
+      aesthetic = "color"
+    ),
+    getAestheticValues(
+      n = totalLength - colorLength,
+      selectionKey = plotConfiguration$points$color,
+      aesthetic = "color"
+    )
+  )
+  
+  # Export the legend captions so the user can update legend keys order
+  plotObject$plotConfiguration$legend$caption <- data.frame(
+    name = colorBreaks,
+    label = colorBreaks,
+    color = colorValues,
+    fill = c(getAestheticValues(n = fillLength, selectionKey = plotConfiguration$ribbons$fill, position = 0, aesthetic = "fill"), 
+             rep(NA, totalLength-fillLength)),
+    linetype = c(getAestheticValues(n = linetypeLength, selectionKey = plotConfiguration$lines$linetype, position = 0, aesthetic = "linetype"), 
+                 rep("blank", totalLength-linetypeLength)),
+    shape = c(rep(" ", totalLength-shapeLength), shapeValues),
+    stringsAsFactors = FALSE
+  )
+  
+  plotObject <- updateTimeProfileLegend(
+    plotObject = plotObject, 
+    caption = plotObject$plotConfiguration$legend$caption
+    ) 
+  
+  if (isIncluded(colorVariable, "legendLabels") & isIncluded(colorObservedVariable, "legendLabels")) {
+    plotObject <- plotObject + ggplot2::guides(color = "none")
   }
+  
+  plotObject <- setLegendPosition(plotObject)
+  plotObject <- setLegendFont(plotObject)
+  eval(parseUpdateAxes())
   return(plotObject)
 }
+
+
+#' @title updateTimeProfileLegend
+#' @description Update time profile legend caption
+#' @param plotObject A ggplot object
+#' @param caption A data.frame as obtained from `getLegendCaption` to use for updating a plot legend.
+#' @return A `ggplot` object
+#' @export
+updateTimeProfileLegend <- function(plotObject, caption){
+  suppressMessages(
+    plotObject <- plotObject + 
+      ggplot2::scale_color_manual(breaks = caption$name, labels = caption$label, values = caption$color) + 
+      ggplot2::guides(
+        fill = "none", shape = "none", linetype = "none",
+        color = ggplot2::guide_legend(
+            override.aes = list(fill = caption$fill, linetype = caption$linetype, shape = caption$shape)
+            )
+          )
+        )
+  plotObject$plotConfiguration$legend$caption <- caption
+  return(plotObject)
+}
\ No newline at end of file
diff --git a/R/timeprofile-datamapping.R b/R/timeprofile-datamapping.R
index dac998d7..f28337ba 100644
--- a/R/timeprofile-datamapping.R
+++ b/R/timeprofile-datamapping.R
@@ -12,6 +12,9 @@ TimeProfileDataMapping <- R6::R6Class(
     #' @param ymin Name of ymin variable to map
     #' @param ymax Name of ymax variable to map
     #' @param group R6 class `Grouping` object or its input
+    #' @param color R6 class `Grouping` object or its input
+    #' @param fill R6 class `Grouping` object or its input
+    #' @param linetype R6 class `Grouping` object or its input
     #' @param data data.frame to map used by `smartMapping`
     #' @return A new `RangeDataMapping` object
     initialize = function(x = NULL,
@@ -19,6 +22,9 @@ TimeProfileDataMapping <- R6::R6Class(
                           ymin = NULL,
                           ymax = NULL,
                           group = NULL,
+                          color = NULL,
+                          fill = NULL,
+                          linetype = NULL,
                           data = NULL) {
 
       # smartMapping is available in utilities-mapping.R
@@ -27,7 +33,7 @@ TimeProfileDataMapping <- R6::R6Class(
         x = x %||% smartMap$x,
         ymin = ymin %||% smartMap$ymin,
         ymax = ymax %||% smartMap$ymax,
-        color = group
+        color = color, fill = fill, linetype = linetype, group = group
       )
       # Since TimeProfileDataMapping inherits from RangeDataMapping
       # super$initialize introduce a self$y which is NULL
@@ -43,7 +49,7 @@ TimeProfileDataMapping <- R6::R6Class(
       validateMapping(self$y, data, nullAllowed = TRUE)
       mapData <- super$checkMapData(data, metaData)
       # This may change depending of how we want to include options
-      if (!isOfLength(self$y, 0)) {
+      if (!isEmpty(self$y)) {
         mapData[, self$y] <- data[, self$y]
       }
       self$data <- mapData
diff --git a/man/ObservedDataMapping.Rd b/man/ObservedDataMapping.Rd
index eaaba3fc..2a42f566 100644
--- a/man/ObservedDataMapping.Rd
+++ b/man/ObservedDataMapping.Rd
@@ -29,8 +29,6 @@ Other DataMapping classes:
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
 \describe{
-\item{\code{lloq}}{mapping lower limit of quantitation variable}
-
 \item{\code{error}}{mapping error bars around scatter points}
 
 \item{\code{mdv}}{mapping missing dependent variable}
@@ -68,8 +66,9 @@ Create a new \code{ObservedDataMapping} object
   error = NULL,
   ymin = NULL,
   ymax = NULL,
-  lloq = NULL,
   mdv = NULL,
+  color = NULL,
+  shape = NULL,
   group = NULL,
   data = NULL
 )}\if{html}{\out{</div>}}
@@ -91,10 +90,12 @@ Deprecated parameter replaced by \code{error}.}
 
 \item{\code{ymax}}{mapping upper end of error bars around scatter points}
 
-\item{\code{lloq}}{mapping lower limit of quantitation variable}
-
 \item{\code{mdv}}{mapping missing dependent variable}
 
+\item{\code{color}}{R6 class \code{Grouping} object or its input}
+
+\item{\code{shape}}{R6 class \code{Grouping} object or its input}
+
 \item{\code{group}}{R6 class \code{Grouping} object or its input}
 
 \item{\code{data}}{data.frame to map used by \code{smartMapping}}
diff --git a/man/RangeDataMapping.Rd b/man/RangeDataMapping.Rd
index b4affc05..615348c1 100644
--- a/man/RangeDataMapping.Rd
+++ b/man/RangeDataMapping.Rd
@@ -65,6 +65,7 @@ Create a new \code{RangeDataMapping} object
   linetype = NULL,
   shape = NULL,
   size = NULL,
+  group = NULL,
   data = NULL
 )}\if{html}{\out{</div>}}
 }
@@ -90,6 +91,8 @@ Create a new \code{RangeDataMapping} object
 
 \item{\code{size}}{R6 class \code{Grouping} object or its input}
 
+\item{\code{group}}{R6 class \code{Grouping} object or its input}
+
 \item{\code{data}}{data.frame to map used by \code{smartMapping}}
 }
 \if{html}{\out{</div>}}
diff --git a/man/TimeProfileDataMapping.Rd b/man/TimeProfileDataMapping.Rd
index c95a379d..83f49dfd 100644
--- a/man/TimeProfileDataMapping.Rd
+++ b/man/TimeProfileDataMapping.Rd
@@ -52,6 +52,9 @@ Create a new \code{TimeProfileDataMapping} object
   ymin = NULL,
   ymax = NULL,
   group = NULL,
+  color = NULL,
+  fill = NULL,
+  linetype = NULL,
   data = NULL
 )}\if{html}{\out{</div>}}
 }
@@ -69,6 +72,12 @@ Create a new \code{TimeProfileDataMapping} object
 
 \item{\code{group}}{R6 class \code{Grouping} object or its input}
 
+\item{\code{color}}{R6 class \code{Grouping} object or its input}
+
+\item{\code{fill}}{R6 class \code{Grouping} object or its input}
+
+\item{\code{linetype}}{R6 class \code{Grouping} object or its input}
+
 \item{\code{data}}{data.frame to map used by \code{smartMapping}}
 }
 \if{html}{\out{</div>}}
diff --git a/man/XYGDataMapping.Rd b/man/XYGDataMapping.Rd
index ef9df8d2..313c192f 100644
--- a/man/XYGDataMapping.Rd
+++ b/man/XYGDataMapping.Rd
@@ -62,6 +62,7 @@ Create a new \code{XYGDataMapping} object
   linetype = NULL,
   shape = NULL,
   size = NULL,
+  group = NULL,
   data = NULL
 )}\if{html}{\out{</div>}}
 }
@@ -85,6 +86,8 @@ Create a new \code{XYGDataMapping} object
 
 \item{\code{size}}{R6 class \code{Grouping} object or its input}
 
+\item{\code{group}}{R6 class \code{Grouping} object or its input}
+
 \item{\code{data}}{data.frame to map used by \code{smartMapping}}
 }
 \if{html}{\out{</div>}}
diff --git a/man/plotTimeProfile.Rd b/man/plotTimeProfile.Rd
index e2c2d19c..0f6ba07d 100644
--- a/man/plotTimeProfile.Rd
+++ b/man/plotTimeProfile.Rd
@@ -25,7 +25,7 @@ plotTimeProfile(
 Unlike \code{data}, meant for simulated data, plotted as lines and ribbons;
 \code{observedData} is plotted as scatter points and errorbars.}
 
-\item{observedDataMapping}{An \code{ObservedDataMapping} object mapping \code{x}, \code{y}, \code{ymin}, \code{ymax}, \code{lloq} and aesthetic groups to their variable names of \code{observedData}.}
+\item{observedDataMapping}{An \code{ObservedDataMapping} object mapping \code{x}, \code{y}, \code{ymin}, \code{ymax} and aesthetic groups to their variable names of \code{observedData}.}
 
 \item{plotConfiguration}{An optional \code{TimeProfilePlotConfiguration} object defining labels, grid, background and watermark.}
 
@@ -42,19 +42,18 @@ Producing Time Profile plots
 obsData <- data.frame(x = c(1, 2, 1, 2, 3), y = c(5, 0.2, 2, 3, 4))
 simTime <- seq(1, 10, 0.1)
 simData <- data.frame(
-x = simTime, 
-y = 10*exp(-simTime),
-ymin = 8*exp(-simTime),
-ymax = 12*exp(-simTime)
+  x = simTime,
+  y = 10 * exp(-simTime),
+  ymin = 8 * exp(-simTime),
+  ymax = 12 * exp(-simTime)
 )
 
 plotTimeProfile(
-data = simData, 
-observedData = obsData,
-dataMapping = TimeProfileDataMapping$new(x = "x", y = "y", ymin = "ymin", ymax = "ymax"),
-observedDataMapping = ObservedDataMapping$new(x = "x", y = "y")
+  data = simData,
+  observedData = obsData,
+  dataMapping = TimeProfileDataMapping$new(x = "x", y = "y", ymin = "ymin", ymax = "ymax"),
+  observedDataMapping = ObservedDataMapping$new(x = "x", y = "y")
 )
-
 }
 \seealso{
 Other molecule plots: 
diff --git a/man/updateTimeProfileLegend.Rd b/man/updateTimeProfileLegend.Rd
new file mode 100644
index 00000000..ba0017a1
--- /dev/null
+++ b/man/updateTimeProfileLegend.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot-timeprofile.R
+\name{updateTimeProfileLegend}
+\alias{updateTimeProfileLegend}
+\title{updateTimeProfileLegend}
+\usage{
+updateTimeProfileLegend(plotObject, caption)
+}
+\arguments{
+\item{plotObject}{A ggplot object}
+
+\item{caption}{A data.frame as obtained from \code{getLegendCaption} to use for updating a plot legend.}
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Update time profile legend caption
+}
diff --git a/vignettes/plot-time-profile.Rmd b/vignettes/plot-time-profile.Rmd
index 080757d5..9559144a 100644
--- a/vignettes/plot-time-profile.Rmd
+++ b/vignettes/plot-time-profile.Rmd
@@ -44,7 +44,7 @@ In such plots, observed data are usually plotted as scatter points with errorbar
 Besides, the usual `tlf` input arguments commonly used by the plot functions (`data`, `metaData`, `dataMapping`, `plotConfiguration` and `plotObject`), the function `plotTimeProfile` also includes the following optional input arguments:
 
 - **`observedData`**: A data.frame to use for plot. Unlike `data`, meant for simulated data, plotted as lines and ribbons; `observedData` is plotted as scatter points and errorbars.
-- **`observedDataMapping`**: An `ObservedDataMapping` object mapping `x`, `y`, `ymin`, `ymax`, `lloq` and aesthetic groups to their variable names of `observedData`.
+- **`observedDataMapping`**: An `ObservedDataMapping` object mapping `x`, `y`, `ymin`, `ymax`, and aesthetic groups to their variable names of `observedData`.
 
 ## 2.2. Data Mappings
 
@@ -76,8 +76,6 @@ As highlighted below, the class is initialized using the method `$new`.
 The optional arguments `ymin` and `ymax` define the error bars of population range or confidence interval.
 The optional argument `error` can be used instead to calculate the errorbars from the variable `y`.
 
-The object also includes the field `lloq` for variables defining the lower limit of quantification of the data. Such data are usually censored and handled differently during modeling.
-
 The argument `mdv`, meant for missing dependent variable, is a flag inspired from Nonmem for excluding data from the final time profile plot.
 
 It can be noted that the argument `group` can be used to group the data by these aesthetics.
@@ -90,7 +88,6 @@ error = NULL,
 ymin = NULL,
 ymax = NULL,
 group = NULL,
-lloq = NULL,
 mdv = NULL
 )
 ```
@@ -401,30 +398,6 @@ plotTimeProfile(
 )
 ```
 
-## 3.2.5. Lower limit of quantification (lloq)
-
-```{r examples lloq}
-# Use common variable before usinf rbind.data.frame
-lloqData <- obsData1
-lloqData$lloq <- 0.5
-
-# Define Data Mapping
-obsDataMapping <- ObservedDataMapping$new(
-  x = "time",
-  y = "concentration",
-  ymin = "minConcentration",
-  ymax = "maxConcentration",
-  group = "caption",
-  lloq = "lloq"
-)
-
-plotTimeProfile(
-  observedData = lloqData,
-  metaData = metaData,
-  observedDataMapping = obsDataMapping
-)
-```
-
 ### 3.2.4. Missing Dependent Variable (mdv)
 
 The following code flags all values higher than `10` as `mdv`, leadin to a plot without any observed data points higher than `10` (removing the first observation)
@@ -547,7 +520,7 @@ obsDataMapping <- ObservedDataMapping$new(
   group = "caption"
 )
 
-plotTimeProfile(
+simAndObsTimeProfile <- plotTimeProfile(
   data = rbind.data.frame(
     simData1,
     simData2
@@ -560,6 +533,43 @@ plotTimeProfile(
   dataMapping = simDataMapping,
   observedDataMapping = obsDataMapping
 )
+simAndObsTimeProfile
+```
+
+In this example observed and simulated data do not share any mapping. 
+Therefore, their color pattern may not match.
+To update the plot and have appropriate colors and labels in the legend, it is possible to get the legend properties and update them as illustrated in the example below that updates the positions of the legend entries:
+
+```{r get legend properties, results='asis'}
+plotLegend <- getLegendCaption(simAndObsTimeProfile) 
+knitr::kable(plotLegend)
+```
+
+```{r update legend properties}
+plotLegend <- plotLegend[c(1,3,2,4), ]
+updateTimeProfileLegend(simAndObsTimeProfile, plotLegend)
+```
+
+Note that if observed and simulated data share common groupings, they will be merged in the final legend.
+
+```{r combine simulated and observed sharing legend}
+commonSimData1 <- simData1
+commonObsData1 <- obsData1
+commonSimData1$caption <- "Common data 1"
+commonObsData1$caption <- "Common data 1"
+commonSimData2 <- simData2
+commonObsData2 <- obsData2
+commonSimData2$caption <- "Common data 2"
+commonObsData2$caption <- "Common data 2"
+
+plotTimeProfile(
+  data = rbind.data.frame(commonSimData1, commonSimData2),
+  observedData = rbind.data.frame(commonObsData1, commonObsData2),
+  metaData = metaData,
+  dataMapping = simDataMapping,
+  observedDataMapping = obsDataMapping
+)
+
 ```
 
 ### 3.3.4. Mutliple simulations and observed data sets with their confidence intervals