diff --git a/DESCRIPTION b/DESCRIPTION
index a7e2deb2..f82aa808 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -100,5 +100,6 @@ Collate:
     'utilities-label.R'
     'utilities-legend.R'
     'utilities-mapping.R'
+    'utilities-molecule-plots.R'
     'utilities-theme.R'
     'utils.R'
diff --git a/R/aaa-utilities.R b/R/aaa-utilities.R
index 69227dfa..187ea083 100644
--- a/R/aaa-utilities.R
+++ b/R/aaa-utilities.R
@@ -42,188 +42,61 @@
   return(parse(text = paste0(objectName, " <- ", value)))
 }
 
-#' @title .parseCheckPlotInputs
-#' @description Create an expression that checks usual plot inputs
-#' @param plotType Type of plot (e.g. "PKRatio" for plotPKRatio)
-#' @return An expression to `eval()`
+#' @title .setDataMapping
+#' @description Set `DataMapping` object internally using `tlf` default if `dataMapping` is not provided
+#' @param dataMapping A `DataMappingClass` object
+#' @param DataMappingClass Required class for `dataMapping`
+#' @param data A data.frame potentially used for smart mapping
+#' @return A `DataMapping` object
 #' @keywords internal
-.parseCheckPlotInputs <- function(plotType) {
-  c(
-    expression(validateIsOfType(data, "data.frame")),
-    expression(validateIsOfType(plotObject, "ggplot", nullAllowed = TRUE)),
-    parse(text = paste0('if(nrow(data)==0){warning(messages$errorNrowData("', plotType, ' plot")); return(plotObject)}')),
-    parse(text = paste0("dataMapping <- dataMapping %||% ", plotType, "DataMapping$new(data=data)")),
-    parse(text = paste0(
-      "plotConfiguration <- plotConfiguration %||% ",
-      plotType, "PlotConfiguration$new(data=data, metaData=metaData, dataMapping=dataMapping)"
-    )),
-    parse(text = paste0('validateIsOfType(dataMapping, "', plotType, 'DataMapping")')),
-    parse(text = paste0('validateIsOfType(plotConfiguration, "', plotType, 'PlotConfiguration")'))
-  )
+.setDataMapping <- function(dataMapping, DataMappingClass, data = NULL) {
+  dataMapping <- dataMapping %||% DataMappingClass$new(data = data)
+  validateIsOfType(dataMapping, DataMappingClass)
+  return(dataMapping)
 }
 
-#' @title .parseUpdateAxes
-#' @description Create an expression that updates the plot axes
-#' @return An expression to `eval()`
+#' @title .setPlotConfiguration
+#' @description Set `PlotConfiguration` object internally using `tlf` default if `plotConfiguration` is not provided
+#' @param plotConfiguration A `PlotConfigurationClass` object
+#' @param PlotConfigurationClass Required class for `plotConfiguration`
+#' @param data A data.frame potentially used for smart plot configuration
+#' @param metaData A list of meta data potentially used for smart plot configuration
+#' @param dataMapping A `DataMapping` object potentially used for smart plot configuration
+#' @return A `PlotConfiguration` object
 #' @keywords internal
-.parseUpdateAxes <- function() {
-  # Try is used to prevent crashes in the final plot due to ggplot2 peculiarities regarding scale functions
-  c(
-    expression(try(suppressMessages(plotObject <- setXAxis(plotObject)))),
-    expression(try(suppressMessages(plotObject <- setYAxis(plotObject))))
-  )
+.setPlotConfiguration <- function(plotConfiguration,
+                                  PlotConfigurationClass,
+                                  data = NULL,
+                                  metaData = NULL,
+                                  dataMapping = NULL) {
+  plotConfiguration <- plotConfiguration %||%
+    PlotConfigurationClass$new(data = data, metaData = metaData, dataMapping = dataMapping)
+  validateIsOfType(plotConfiguration, PlotConfigurationClass)
+  return(plotConfiguration)
 }
 
-#' @title .parseUpdateAestheticProperty
-#' @description Create an expression that updates the aesthetic properties based on
-#' the information of `PlotConfiguration`
-#' @param aestheticProperty Name of aesthetic property as defined in `AestheticProperties`
-#' @param plotConfigurationProperty Name of PlotConfiguration property as defined in `AestheticProperties`
-#' @return An expression to `eval()`
+#' @title .setPlotObject
+#' @description Set a `ggplot` object associated with its `plotConfiguration`
+#' @param plotObject A `ggplot` object
+#' @param plotConfiguration A `PlotConfiguration` object
+#' @return A `ggplot` object
 #' @keywords internal
-.parseUpdateAestheticProperty <- function(aestheticProperty, plotConfigurationProperty) {
-  c(
-    parse(text = paste0(aestheticProperty, 'Variable <- gsub("`", "", mapLabels$', aestheticProperty, ")")),
-    parse(text = paste0(aestheticProperty, "Length <- length(unique(mapData[, ", aestheticProperty, "Variable]))")),
-    # Update the property using ggplot `scale` functions
-    parse(text = paste0(
-      "suppressMessages(plotObject <- plotObject + ggplot2::scale_", aestheticProperty, "_manual(",
-      "values=.getAestheticValues(n=", aestheticProperty, "Length,",
-      "selectionKey=plotConfiguration$", plotConfigurationProperty, "$", 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')}"))
-  )
+.setPlotObject <- function(plotObject,
+                           plotConfiguration = NULL) {
+  plotObject <- plotObject %||% initializePlot(plotConfiguration)
+  validateIsOfType(plotObject, "ggplot")
+  validateIsIncluded("plotConfiguration", names(plotObject))
+  return(plotObject)
 }
 
-#' @title .parseAddScatterLayer
-#' @description Create an expression that adds scatter plot layer
-#' @return An expression to `eval()`
-#' @keywords internal
-.parseAddScatterLayer <- function() {
-  expression({
-    plotObject <- plotObject +
-      ggplot2::geom_point(
-        data = mapData,
-        mapping = ggplot2::aes_string(
-          x = mapLabels$x,
-          y = mapLabels$y,
-          color = mapLabels$color,
-          shape = mapLabels$shape
-        ),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$size, position = 0, aesthetic = "size"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$alpha, aesthetic = "alpha"),
-        na.rm = TRUE,
-        show.legend = TRUE
-      )
-  })
-}
 
-#' @title .parseAddLineLayer
-#' @description Create an expression that adds scatter plot layer
-#' TODO: create a vignette explaining how argument `lines` in dataMapping is related to this
-#' @param type one of "horizontal", "vertical" or "diagonal"
-#' Note that for "diagonal", geom_abline is used.
-#' `value` of intercept is taken as is for linear scale but corresponds to the log of `value` for log scale.
-#'  For instance, intercept = c(-1, 0, 1) with log scale actually means that the line will go through c(0.1, 1, 10)
-#'  because c(-1, 0, 1) = log10(c(0.1, 1, 10)).
-#' @param value value of xintercept or yintercept
-#' @param position line position for aesthetic properties
-#' @return An expression to `eval()`
-#' @keywords internal
-.parseAddLineLayer <- function(type, value, position) {
-  parse(text = paste0(
-    "plotObject <- plotObject + ",
-    switch(type,
-      "horizontal" = paste0("ggplot2::geom_hline(yintercept = ", value, ","),
-      "vertical" = paste0("ggplot2::geom_vline(xintercept = ", value, ","),
-      "diagonal" = paste0("ggplot2::geom_abline(slope=1, intercept = ", value, ","),
-      "ddiHorizontal" = paste0("ggplot2::geom_abline(slope=0, intercept = ", value, ",")
-    ),
-    "color=.getAestheticValues(n=1,selectionKey=plotConfiguration$lines$color,position=", position, ',aesthetic="color"),',
-    "linetype=.getAestheticValues(n=1,selectionKey=plotConfiguration$lines$linetype,position=", position, ',aesthetic="linetype"),',
-    "alpha=.getAestheticValues(n=1,selectionKey=plotConfiguration$lines$alpha,position=", position, ',aesthetic="alpha"),',
-    "size=.getAestheticValues(n=1,selectionKey=plotConfiguration$lines$size,position=", position, ', aesthetic="size"))'
-  ))
-}
-
-#' @title .parseAddUncertaintyLayer
-#' @description Create an expression that adds errorbars
-#' `mapLabels` needs to be obtained from `DataMapping` objects
-#' @return An expression to `eval()`
+#' @title .updateAxes
+#' @description Updates the plot axes
+#' @return A `ggplot` object
 #' @keywords internal
-.parseAddUncertaintyLayer <- function(direction = "vertical") {
-  parse(text = paste0(
-    "plotObject <- plotObject +",
-    # Plot error bars from xmin/ymin to x/y
-    # If lower value is negative and plot is log scaled,
-    # Upper bar will still be plotted
-    "ggplot2::geom_linerange(",
-    "data = mapData,",
-    "mapping = aes_string(",
-    switch(direction,
-      "vertical" = "x = mapLabels$x, ymin = mapLabels$ymin, ymax = mapLabels$y,",
-      "horizontal" = "y = mapLabels$y, xmin = mapLabels$xmin, xmax = mapLabels$x,"
-    ),
-    "color = mapLabels$color,",
-    "group = mapLabels$color",
-    "),",
-    'size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$size, position = 0, aesthetic = "size"),',
-    'linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$linetype, aesthetic = "linetype"),',
-    'alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, aesthetic = "alpha"),',
-    "na.rm = TRUE,",
-    "show.legend = FALSE",
-    ") + ",
-    "ggplot2::geom_linerange(",
-    "data = mapData,",
-    "mapping = aes_string(",
-    switch(direction,
-      "vertical" = "x = mapLabels$x, ymin = mapLabels$y, ymax = mapLabels$ymax,",
-      "horizontal" = "y = mapLabels$y, xmin = mapLabels$x, xmax = mapLabels$xmax,"
-    ),
-    "color = mapLabels$color,",
-    "group = mapLabels$color",
-    "),",
-    'size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$size, position = 0, aesthetic = "size"),',
-    'linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$linetype, aesthetic = "linetype"),',
-    'alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, aesthetic = "alpha"),',
-    "na.rm = TRUE,",
-    "show.legend = FALSE",
-    ") + ",
-    # Add lower cap to error bar
-    "ggplot2::geom_point(",
-    "data = mapData,",
-    "mapping = aes_string(",
-    switch(direction,
-           "vertical" = "x = mapLabels$x, y = mapLabels$ymin,",
-           "horizontal" = "y = mapLabels$y, x = mapLabels$xmin,"
-    ),
-    "color = mapLabels$color,",
-    "group = mapLabels$color",
-    "),",
-    'size = tlfEnv$defaultErrorbarCapSize,',
-    'shape = ', switch(direction, "vertical" = '"_"', "horizontal" = '"|"'), ',',
-    'alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, aesthetic = "alpha"),',
-    "na.rm = TRUE,",
-    "show.legend = FALSE",
-    ") + ",
-    # Add upper cap to error bar
-    "ggplot2::geom_point(",
-    "data = mapData,",
-    "mapping = aes_string(",
-    switch(direction,
-           "vertical" = "x = mapLabels$x, y = mapLabels$ymax,",
-           "horizontal" = "y = mapLabels$y, x = mapLabels$xmax,"
-    ),
-    "color = mapLabels$color,",
-    "group = mapLabels$color",
-    "),",
-    'size = tlfEnv$defaultErrorbarCapSize,',
-    'shape = ', switch(direction, "vertical" = '"_"', "horizontal" = '"|"'), ',', 
-    'alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, aesthetic = "alpha"),',
-    "na.rm = TRUE,",
-    "show.legend = FALSE",
-    ")"
-  ))
+.updateAxes <- function(plotObject) {
+  # Try is used to prevent crashes in the final plot due to ggplot2 peculiarities regarding scale functions
+  try(suppressMessages(plotObject <- setXAxis(plotObject)))
+  try(suppressMessages(plotObject <- setYAxis(plotObject)))
+  return(plotObject)
 }
diff --git a/R/atom-plots.R b/R/atom-plots.R
index e1209cdb..884f98e3 100644
--- a/R/atom-plots.R
+++ b/R/atom-plots.R
@@ -184,7 +184,7 @@ addScatter <- function(data = NULL,
     newLabels = newLabels,
     aestheticSelections = plotConfiguration$points
   ))
-  eval(.parseUpdateAxes())
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
@@ -362,7 +362,7 @@ addLine <- function(data = NULL,
     newLabels = newLabels,
     aestheticSelections = plotConfiguration$lines
   ))
-  eval(.parseUpdateAxes())
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
@@ -495,7 +495,7 @@ addRibbon <- function(data = NULL,
     newLabels = newLabels,
     aestheticSelections = plotConfiguration$ribbons
   ))
-  eval(.parseUpdateAxes())
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
@@ -566,7 +566,7 @@ addErrorbar <- function(data = NULL,
                         plotObject = NULL) {
   validateIsOfType(dataMapping, c("RangeDataMapping", "ObservedDataMapping"), nullAllowed = TRUE)
   validateIsOfType(plotConfiguration, PlotConfiguration, nullAllowed = TRUE)
-  
+
   # If data is not input, creates data and its mapping from x, ymin and ymax input
   if (isEmpty(data)) {
     data <- as.data.frame(cbind(x = x, ymin = ymin %||% 0, ymax = ymax %||% 0))
@@ -600,31 +600,31 @@ addErrorbar <- function(data = NULL,
   }
   mapData$legendLabels <- caption %||% mapData$legendLabels
   legendLength <- length(unique(mapData$legendLabels))
-  
+
   eval(.parseVariableToObject("plotObject$plotConfiguration$errorbars", c("color", "size", "linetype"), keepIfNull = TRUE))
-  
+
   plotObject <- plotObject +
     ggplot2::geom_linerange(
       data = mapData,
       mapping = aes_string(
-        x = mapLabels$x, 
-        ymin = mapLabels$ymin, 
+        x = mapLabels$x,
+        ymin = mapLabels$ymin,
         ymax = mapLabels$ymax,
         color = mapLabels$color,
         group = mapLabels$color
-        ),
+      ),
       size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$size, position = 0, aesthetic = "size"),
       linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$linetype, aesthetic = "linetype"),
       alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, aesthetic = "alpha"),
       color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$color, aesthetic = "color"),
       na.rm = TRUE,
       show.legend = FALSE
-      ) + 
+    ) +
     ggplot2::geom_point(
       data = mapData,
       mapping = aes_string(
-        x = mapLabels$x, 
-        y = mapLabels$ymin, 
+        x = mapLabels$x,
+        y = mapLabels$ymin,
         color = mapLabels$color,
         group = mapLabels$color
       ),
@@ -634,11 +634,11 @@ addErrorbar <- function(data = NULL,
       color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$color, aesthetic = "color"),
       na.rm = TRUE,
       show.legend = FALSE
-    ) + 
+    ) +
     ggplot2::geom_point(
       data = mapData,
       mapping = aes_string(
-        x = mapLabels$x, 
+        x = mapLabels$x,
         y = mapLabels$ymax,
         color = mapLabels$color,
         group = mapLabels$color
@@ -650,9 +650,9 @@ addErrorbar <- function(data = NULL,
       na.rm = TRUE,
       show.legend = FALSE
     )
-  
+
   # Try is used to prevent crashes in the final plot due to ggplot2 peculiarities regarding scale functions
-  eval(.parseUpdateAxes())
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
diff --git a/R/datamapping-range.R b/R/datamapping-range.R
index 1c4a1884..34050241 100644
--- a/R/datamapping-range.R
+++ b/R/datamapping-range.R
@@ -67,26 +67,34 @@ RangeDataMapping <- R6::R6Class(
 
       # All possible Groupings are listed in the enum LegendTypes
       for (groupType in LegendTypes) {
-        if (!is.null(self$groupMapping[[groupType]]$group)) {
-          grouping <- self$groupMapping[[groupType]]
+        if (isEmpty(self$groupMapping[[groupType]]$group)) {
+          next
+        }
+        grouping <- self$groupMapping[[groupType]]
+
+        groupVariables <- grouping$group
+        if (isOfType(groupVariables, "data.frame")) {
+          # Last group variable is the label in group data.frames
+          # and need to be removed from the check
+          groupVariables <- names(groupVariables)
+          groupVariables <- utils::head(groupVariables, -1)
+        }
+        .validateMapping(groupVariables, data)
+        # Enforce grouping variables to be factors
+        self$data[, grouping$label] <- as.factor(grouping$getCaptions(data, metaData))
+
+        # Dummy variable for default aesthetics that will be used to define legend labels
+        legendLabels <- self$data$legendLabels %||% grouping$getCaptions(data, metaData)
 
-          groupVariables <- grouping$group
-          if (isOfType(groupVariables, "data.frame")) {
-            # Last group variable is the label in group data.frames
-            # and need to be removed from the check
-            groupVariables <- names(groupVariables)
-            groupVariables <- utils::head(groupVariables, -1)
-          }
-          .validateMapping(groupVariables, data)
-          self$data[, grouping$label] <- grouping$getCaptions(data, metaData)
-          # Dummy variable for default aesthetics
-          # Will be used to define legend labels
-          self$data$legendLabels <- ifNotNull(
-            self$data$legendLabels,
-            paste(self$data$legendLabels, grouping$getCaptions(data, metaData), sep = "-"),
-            grouping$getCaptions(data, metaData)
-          )
+        # Prevent duplication of legend if groupings are the same
+        if (isTRUE(all.equal(legendLabels, grouping$getCaptions(data, metaData)))) {
+          self$data$legendLabels <- legendLabels
+          next
         }
+        self$data$legendLabels <- as.factor(paste(as.character(self$data$legendLabels),
+          as.character(grouping$getCaptions(data, metaData)),
+          sep = "-"
+        ))
       }
 
       if (is.null(self$data$legendLabels)) {
diff --git a/R/datamapping-xygroup.R b/R/datamapping-xygroup.R
index a456f588..2396cc64 100644
--- a/R/datamapping-xygroup.R
+++ b/R/datamapping-xygroup.R
@@ -68,7 +68,7 @@ XYGDataMapping <- R6::R6Class(
 
       # All possible Groupings are listed in the enum LegendTypes
       for (groupType in LegendTypes) {
-        if (isOfLength(self$groupMapping[[groupType]]$group, 0)) {
+        if (isEmpty(self$groupMapping[[groupType]]$group)) {
           next
         }
         grouping <- self$groupMapping[[groupType]]
diff --git a/R/plot-boxwhisker.R b/R/plot-boxwhisker.R
index 577d7905..751e8b97 100644
--- a/R/plot-boxwhisker.R
+++ b/R/plot-boxwhisker.R
@@ -34,36 +34,25 @@ plotBoxWhisker <- function(data,
                            dataMapping = NULL,
                            plotConfiguration = NULL,
                            plotObject = NULL) {
-  dataMapping <- dataMapping %||% BoxWhiskerDataMapping$new(data = data)
-  plotConfiguration <- plotConfiguration %||% BoxWhiskerPlotConfiguration$new(
-    data = data,
-    metaData = metaData,
-    dataMapping = dataMapping
-  )
-
-  validateIsOfType(dataMapping, "BoxWhiskerDataMapping")
-  validateIsOfType(plotConfiguration, "BoxWhiskerPlotConfiguration")
+  #----- Validation and formatting of input arguments -----
+  validateIsNotEmpty(data)
   validateIsOfType(data, "data.frame")
-  validateIsOfType(plotObject, "ggplot", nullAllowed = TRUE)
-
+  dataMapping <- .setDataMapping(dataMapping, BoxWhiskerDataMapping, data)
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, BoxWhiskerPlotConfiguration,
+    data, metaData, dataMapping
+  )
   # Overwrites plotConfiguration$outliers if outliers is not null
   validateIsLogical(outliers, nullAllowed = TRUE)
   plotConfiguration$outliers <- outliers
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
 
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
-  if (nrow(data) == 0) {
-    warning(messages$errorNrowData("box whisker plot"))
-    return(plotObject)
-  }
-
-  # Add Plot Configuration layers and box whisker plots
+  #----- Build layers of molecule plot -----
   plotObject <- .addBoxWhisker(data, metaData, dataMapping, plotConfiguration, plotObject)
   if (plotConfiguration$outliers) {
     plotObject <- .addOutliers(data, metaData, dataMapping, plotConfiguration, plotObject)
   }
-  try(suppressMessages(plotObject <- setXAxis(plotObject)))
-  try(suppressMessages(plotObject <- setYAxis(plotObject)))
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
@@ -80,6 +69,13 @@ plotBoxWhisker <- function(data,
   # Convert the mapping into characters usable by aes_string
   mapLabels <- .getAesStringMapping(dataMapping)
 
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = 0,
+    plotConfigurationProperty = plotObject$plotConfiguration$ribbons,
+    propertyNames = c("size", "alpha", "linetype")
+  )
+
   plotObject <- plotObject +
     ggplot2::geom_boxplot(
       data = mapData,
@@ -93,30 +89,20 @@ plotBoxWhisker <- function(data,
         fill = mapLabels$fill,
         color = mapLabels$color
       ),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$alpha, position = 0, aesthetic = "alpha"),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$size, position = 0, aesthetic = "size"),
-      linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$linetype, position = 0, aesthetic = "linetype"),
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size,
+      linetype = aestheticValues$linetype,
       show.legend = TRUE,
       stat = "identity"
     )
 
-  # Define linetype, color, f# Define shapes and colors based on plotConfiguration$points properties
-  fillVariable <- gsub("`", "", mapLabels$fill)
-  colorVariable <- gsub("`", "", mapLabels$color)
-  fillLength <- length(unique(mapData[, fillVariable]))
-  colorLength <- length(unique(mapData[, colorVariable]))
-
-  plotObject <- plotObject +
-    ggplot2::scale_fill_manual(values = .getAestheticValues(n = fillLength, selectionKey = plotConfiguration$ribbons$fill, aesthetic = "fill")) +
-    ggplot2::scale_color_manual(values = .getAestheticValues(n = colorLength, selectionKey = plotConfiguration$ribbons$color, aesthetic = "color"))
-
-  # If variable is legendLabel, remove it from legend
-  if (isIncluded(fillVariable, "legendLabels")) {
-    plotObject <- plotObject + ggplot2::guides(fill = "none")
-  }
-  if (isIncluded(colorVariable, "legendLabels")) {
-    plotObject <- plotObject + ggplot2::guides(color = "none")
-  }
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "ribbons",
+    propertyNames = c("color", "fill"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
   return(plotObject)
 }
 
@@ -132,7 +118,14 @@ plotBoxWhisker <- function(data,
   # Convert the mapping into characters usable by aes_string
   mapLabels <- .getAesStringMapping(dataMapping)
 
-  # addScatter cannot be used in this case,
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = 0,
+    plotConfigurationProperty = plotObject$plotConfiguration$points,
+    propertyNames = c("size", "alpha", "shape", "color")
+  )
+
+  # addScatterLayer cannot be used in this case,
   # because position dodge is needed to align boxes and outlier points
   # no matter the number of groups, the value of 0.9 will be always fix
   # otherwise, points won't be centered anymore
@@ -147,9 +140,10 @@ plotBoxWhisker <- function(data,
         group = mapLabels$fill,
         color = mapLabels$color
       ),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$size, position = 0, aesthetic = "size"),
-      shape = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$shape, position = 0, aesthetic = "shape"),
-      color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$color, position = 0, aesthetic = "color"),
+      size = aestheticValues$size,
+      shape = aestheticValues$shape,
+      color = aestheticValues$color,
+      alpha = aestheticValues$alpha,
       show.legend = TRUE,
       na.rm = TRUE,
       position = position_dodge(width = 0.9)
@@ -162,9 +156,10 @@ plotBoxWhisker <- function(data,
         group = mapLabels$fill,
         color = mapLabels$color
       ),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$size, position = 0, aesthetic = "size"),
-      shape = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$shape, position = 0, aesthetic = "shape"),
-      color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$color, position = 0, aesthetic = "color"),
+      size = aestheticValues$size,
+      shape = aestheticValues$shape,
+      color = aestheticValues$color,
+      alpha = aestheticValues$alpha,
       show.legend = TRUE,
       na.rm = TRUE,
       position = position_dodge(width = 0.9)
diff --git a/R/plot-ddiratio.R b/R/plot-ddiratio.R
index 35013465..f599e373 100644
--- a/R/plot-ddiratio.R
+++ b/R/plot-ddiratio.R
@@ -44,64 +44,91 @@ plotDDIRatio <- function(data,
                          foldDistance = NULL,
                          deltaGuest = NULL,
                          plotObject = NULL) {
-  eval(.parseCheckPlotInputs("DDIRatio"))
-  validateIsLogical(residualsVsObserved, nullAllowed = TRUE)
+  #----- Validation and formatting of input arguments -----
+  validateIsNotEmpty(data)
+  validateIsOfType(data, "data.frame")
+  dataMapping <- .setDataMapping(dataMapping, DDIRatioDataMapping, data)
   validateIsNumeric(foldDistance, nullAllowed = TRUE)
   validateIsNumeric(deltaGuest, nullAllowed = TRUE)
-  mapData <- dataMapping$checkMapData(data)
-  mapLabels <- .getAesStringMapping(dataMapping)
-
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
+  validateIsLogical(residualsVsObserved, nullAllowed = TRUE)
+  if (!isEmpty(foldDistance)) {
+    dataMapping$lines <- getLinesFromFoldDistance(foldDistance)
+  }
   dataMapping$residualsVsObserved <- residualsVsObserved %||% dataMapping$residualsVsObserved
   dataMapping$deltaGuest <- deltaGuest %||% dataMapping$deltaGuest
 
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, DDIRatioPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
+
+  mapData <- dataMapping$checkMapData(data)
+  mapLabels <- .getAesStringMapping(dataMapping)
+
+  #----- Build layers of molecule plot -----
+  # Each new layer is added on top of previous
+  # Thus, scatter points are added as last layer to prevent them being hidden by lines or errorbars
+  # 1- Horizontal/Diagonal lines
   lineOrientation <- "diagonal"
   if (dataMapping$residualsVsObserved) {
     lineOrientation <- "ddiHorizontal"
   }
-  # Include diagonal or horizontal lines depending on the plot type
-  if (!isEmpty(foldDistance)) {
-    dataMapping$lines <- getLinesFromFoldDistance(foldDistance)
-  }
   for (lineIndex in seq_along(dataMapping$lines)) {
     lineValue <- .getAblineValues(dataMapping$lines[[lineIndex]], plotConfiguration$yAxis$scale)
-    # position correspond to the number of layer lines already added
-    eval(.parseAddLineLayer(lineOrientation, lineValue, lineIndex - 1))
-  }
-  if (isEmpty(lineIndex)) {
-    lineIndex <- 0
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = lineOrientation,
+      value = lineValue,
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
-  # Add Guest et al. lines to plot
+  lineIndex <- ifNotNull(lineIndex, lineIndex, 0)
+
+  # 2- Guest et al. lines
   # guestData is a data.frame with x, ymin and ymax
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = lineIndex,
+    plotConfigurationProperty = plotObject$plotConfiguration$lines,
+    propertyNames = c("color", "linetype", "size", "alpha")
+  )
   guestData <- getGuestValuesFromDataMapping(data, dataMapping)
   plotObject <- plotObject +
     ggplot2::geom_path(
       data = guestData,
       mapping = ggplot2::aes_string(x = "x", y = "ymin"),
       na.rm = TRUE,
-      color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-      linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
     ) +
     ggplot2::geom_path(
       data = guestData,
       mapping = ggplot2::aes_string(x = "x", y = "ymax"),
       na.rm = TRUE,
-      color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-      linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
     )
 
+  # 2- Error bars
+  plotObject <- .addErrorbarLayer(plotObject, data = mapData, mapLabels = mapLabels)
+  # 3- Scatter points
+  plotObject <- .addScatterLayer(plotObject, data = mapData, mapLabels = mapLabels)
 
-  eval(.parseAddUncertaintyLayer())
-  eval(.parseAddScatterLayer())
-  # Define shapes and colors based on plotConfiguration$points properties
-  eval(.parseUpdateAestheticProperty(AestheticProperties$color, "points"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$shape, "points"))
-  eval(.parseUpdateAxes())
+  #----- Update properties using ggplot2::scale functions -----
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "points",
+    propertyNames = c("color", "shape"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
diff --git a/R/plot-histogram.R b/R/plot-histogram.R
index 43f4b3d3..36df0cf1 100644
--- a/R/plot-histogram.R
+++ b/R/plot-histogram.R
@@ -43,11 +43,7 @@ plotHistogram <- function(data = NULL,
                           distribution = NULL,
                           plotConfiguration = NULL,
                           plotObject = NULL) {
-  validateIsNumeric(bins, nullAllowed = TRUE)
-  validateIsNumeric(binwidth, nullAllowed = TRUE)
-  validateIsLogical(stack, nullAllowed = TRUE)
-  validateIsIncluded(distribution, c("normal", "logNormal", "none"), nullAllowed = TRUE)
-
+  #----- Validation and formatting of input arguments -----
   if (is.null(data)) {
     validateIsNumeric(x)
     data <- data.frame(x = x)
@@ -56,45 +52,50 @@ plotHistogram <- function(data = NULL,
       data = data
     )
   }
-  dataMapping <- dataMapping %||% HistogramDataMapping$new(data = data)
-  plotConfiguration <- plotConfiguration %||% HistogramPlotConfiguration$new(
-    data = data,
-    metaData = metaData,
-    dataMapping = dataMapping
-  )
-
-  validateIsOfType(dataMapping, "HistogramDataMapping")
-  validateIsOfType(plotConfiguration, "HistogramPlotConfiguration")
+  validateIsNotEmpty(data)
   validateIsOfType(data, "data.frame")
-  validateIsOfType(plotObject, "ggplot", nullAllowed = TRUE)
+  dataMapping <- .setDataMapping(dataMapping, HistogramDataMapping, data)
+
+  # Update dataMapping if inputs provided by user
+  validateIsNumeric(bins, nullAllowed = TRUE)
+  validateIsNumeric(binwidth, nullAllowed = TRUE)
+  validateIsLogical(stack, nullAllowed = TRUE)
+  validateIsIncluded(distribution, c("normal", "logNormal", "none"), nullAllowed = TRUE)
 
-  # Overwrites plotConfiguration and dataMapping if some inputs are not null
   dataMapping$stack <- stack %||% dataMapping$stack
   dataMapping$distribution <- distribution %||% dataMapping$distribution
   dataMapping$bins <- bins %||% dataMapping$bins
   dataMapping$binwidth <- binwidth %||% dataMapping$binwidth
 
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
-  if (nrow(data) == 0) {
-    warning(messages$errorNrowData("Histogram"))
-    return(plotObject)
-  }
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, HistogramPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
 
-  # Get transformed data from mapping and convert labels into characters usable by aes_string
   mapData <- dataMapping$checkMapData(data)
   mapLabels <- .getAesStringMapping(dataMapping)
 
+  #----- Build layers of molecule plot -----
+  # position defines if bars are stacked or plotted side by side
   position <- ggplot2::position_nudge()
   if (dataMapping$stack) {
     position <- ggplot2::position_stack()
   }
 
+  # If argument bins is of length > 1,
+  # bins corresponds to bin edges instead of number of bins
   edges <- NULL
   if (length(dataMapping$bins) > 1) {
     edges <- dataMapping$bins
   }
 
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    plotConfigurationProperty = plotObject$plotConfiguration$ribbons,
+    propertyNames = c("color", "size", "alpha", "linetype")
+  )
+  # 1- Histogram
   plotObject <- plotObject +
     ggplot2::geom_histogram(
       data = mapData,
@@ -106,16 +107,24 @@ plotHistogram <- function(data = NULL,
       bins = dataMapping$bins,
       binwidth = dataMapping$binwidth,
       breaks = edges,
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$size, position = 0, aesthetic = "size"),
-      color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$color, position = 0, aesthetic = "color"),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$alpha, position = 0, aesthetic = "alpha")
+      size = aestheticValues$size,
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha
     )
 
   # If distribution is provided by dataMapping, get median and distribution of the data
   fitData <- .getDistributionFit(mapData, dataMapping)
   fitMedian <- .getDistributionMed(mapData, dataMapping)
 
-  if (!isOfLength(fitData, 0)) {
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    plotConfigurationProperty = plotObject$plotConfiguration$lines,
+    propertyNames = c("size", "alpha")
+  )
+
+  # 2- Lines of distribution fit
+  if (!isEmpty(fitData)) {
     plotObject <- plotObject +
       ggplot2::geom_line(
         data = fitData,
@@ -125,21 +134,38 @@ plotHistogram <- function(data = NULL,
           color = "legendLabels",
           linetype = "legendLabels"
         ),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = 0, aesthetic = "size")
+        size = aestheticValues$size,
+        alpha = aestheticValues$alpha
       )
   }
 
-  # Include vertical lines
+  # 3- Vertical lines of median
   for (lineIndex in seq_along(fitMedian)) {
-    # position corresponds to the number of layer lines already added
-    eval(.parseAddLineLayer("vertical", fitMedian[lineIndex], lineIndex - 1))
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = "vertical",
+      value = fitMedian[lineIndex],
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
 
-  # Define fill based on plotConfiguration$points properties
-  eval(.parseUpdateAestheticProperty(AestheticProperties$fill, "ribbons"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$color, "lines"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$linetype, "lines"))
-  eval(.parseUpdateAxes())
+  #----- Update properties using ggplot2::scale functions -----
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "ribbons",
+    propertyNames = "fill",
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "lines",
+    propertyNames = c("color", "linetype"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
diff --git a/R/plot-obs-vs-pred.R b/R/plot-obs-vs-pred.R
index 985e024e..e4329a52 100644
--- a/R/plot-obs-vs-pred.R
+++ b/R/plot-obs-vs-pred.R
@@ -51,28 +51,49 @@ plotObsVsPred <- function(data,
                           foldDistance = NULL,
                           smoother = NULL,
                           plotObject = NULL) {
-  eval(.parseCheckPlotInputs("ObsVsPred"))
-  validateIsIncluded(smoother, c("loess", "lm"), nullAllowed = TRUE)
+  #----- Validation and formatting of input arguments -----
+  validateIsNotEmpty(data)
+  validateIsOfType(data, "data.frame")
+  dataMapping <- .setDataMapping(dataMapping, ObsVsPredDataMapping, data)
   validateIsNumeric(foldDistance, nullAllowed = TRUE)
+  if (!isEmpty(foldDistance)) {
+    dataMapping$lines <- getLinesFromFoldDistance(foldDistance)
+  }
+  validateIsIncluded(smoother, c("loess", "lm"), nullAllowed = TRUE)
   dataMapping$smoother <- smoother %||% dataMapping$smoother
+
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, ObsVsPredPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
+
   mapData <- dataMapping$checkMapData(data)
   mapLabels <- .getAesStringMapping(dataMapping)
 
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
-  # Add diagonal lines with offset defined in lines of dataMapping
-  if (!isEmpty(foldDistance)) {
-    dataMapping$lines <- getLinesFromFoldDistance(foldDistance)
-  }
+  #----- Build layers of molecule plot -----
+  # Each new layer is added on top of previous
+  # Thus, scatter points are added as last layer to prevent them being hidden by lines or errorbars
+  # 1- Diagonal lines
   for (lineIndex in seq_along(dataMapping$lines)) {
     lineValue <- .getAblineValues(dataMapping$lines[[lineIndex]], plotConfiguration$yAxis$scale)
-    # position correspond to the number of layer lines already added
-    eval(.parseAddLineLayer("diagonal", lineValue, lineIndex - 1))
-  }
-  if (isEmpty(lineIndex)) {
-    lineIndex <- 0
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = "diagonal",
+      value = lineValue,
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
-  # Add Smoother if defined
+  lineIndex <- ifNotNull(lineIndex, lineIndex, 0)
+
+  # 2- Smoother line
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = lineIndex,
+    plotConfigurationProperty = plotObject$plotConfiguration$lines,
+    propertyNames = c("color", "linetype", "size", "alpha")
+  )
   if (isIncluded(dataMapping$smoother, "loess")) {
     plotObject <- plotObject +
       ggplot2::geom_smooth(
@@ -82,10 +103,10 @@ plotObsVsPred <- function(data,
         formula = "y ~ x",
         na.rm = TRUE,
         mapping = ggplot2::aes_string(x = mapLabels$x, y = mapLabels$y),
-        color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-        linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+        color = aestheticValues$color,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        size = aestheticValues$size
       )
   }
   if (isIncluded(dataMapping$smoother, "lm")) {
@@ -97,18 +118,25 @@ plotObsVsPred <- function(data,
         se = FALSE,
         formula = "y ~ x",
         na.rm = TRUE,
-        color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-        linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+        color = aestheticValues$color,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        size = aestheticValues$size
       )
   }
+  # 2- Error bars
+  plotObject <- .addErrorbarLayer(plotObject, data = mapData, mapLabels = mapLabels, direction = "horizontal")
+  # 3- Scatter points
+  plotObject <- .addScatterLayer(plotObject, data = mapData, mapLabels = mapLabels)
 
-  eval(.parseAddUncertaintyLayer(direction = "horizontal"))
-  eval(.parseAddScatterLayer())
-  # Define shapes and colors based on plotConfiguration$points properties
-  eval(.parseUpdateAestheticProperty(AestheticProperties$color, "points"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$shape, "points"))
-  eval(.parseUpdateAxes())
+  #----- Update properties using ggplot2::scale functions -----
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "points",
+    propertyNames = c("color", "shape"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
diff --git a/R/plot-pkratio.R b/R/plot-pkratio.R
index 4da152cf..0eae1e43 100644
--- a/R/plot-pkratio.R
+++ b/R/plot-pkratio.R
@@ -37,28 +37,49 @@ plotPKRatio <- function(data,
                         plotConfiguration = NULL,
                         foldDistance = NULL,
                         plotObject = NULL) {
-  eval(.parseCheckPlotInputs("PKRatio"))
+  #----- Validation and formatting of input arguments -----
+  validateIsNotEmpty(data)
+  validateIsOfType(data, "data.frame")
+  dataMapping <- .setDataMapping(dataMapping, PKRatioDataMapping, data)
   validateIsNumeric(foldDistance, nullAllowed = TRUE)
-  mapData <- dataMapping$checkMapData(data)
-  mapLabels <- .getAesStringMapping(dataMapping)
-
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
-  # Include horizontal lines
   if (!isEmpty(foldDistance)) {
     dataMapping$lines <- getLinesFromFoldDistance(foldDistance)
   }
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, PKRatioPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
+
+  mapData <- dataMapping$checkMapData(data)
+  mapLabels <- .getAesStringMapping(dataMapping)
+
+  #----- Build layers of molecule plot -----
+  # Each new layer is added on top of previous
+  # Thus, scatter points are added as last layer to prevent them being hidden by lines or errorbars
+  # 1- Horizontal lines
   for (lineIndex in seq_along(dataMapping$lines)) {
-    # position correspond to the number of layer lines already added
-    eval(.parseAddLineLayer("horizontal", dataMapping$lines[[lineIndex]], lineIndex - 1))
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = "horizontal",
+      value = dataMapping$lines[[lineIndex]],
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
+  # 2- Error bars
+  plotObject <- .addErrorbarLayer(plotObject, data = mapData, mapLabels = mapLabels)
+  # 3- Scatter points
+  plotObject <- .addScatterLayer(plotObject, data = mapData, mapLabels = mapLabels)
 
-  # If uncertainty is defined, add error bars
-  eval(.parseAddUncertaintyLayer())
-  eval(.parseAddScatterLayer())
-  # Define shapes and colors based on plotConfiguration$points properties
-  eval(.parseUpdateAestheticProperty(AestheticProperties$color, "points"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$shape, "points"))
-  eval(.parseUpdateAxes())
+  #----- Update properties using ggplot2::scale functions -----
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "points",
+    propertyNames = c("color", "shape"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
diff --git a/R/plot-res-vs-pred.R b/R/plot-res-vs-pred.R
index 786386c6..2d285212 100644
--- a/R/plot-res-vs-pred.R
+++ b/R/plot-res-vs-pred.R
@@ -29,22 +29,44 @@ plotResVsPred <- function(data,
                           plotConfiguration = NULL,
                           smoother = NULL,
                           plotObject = NULL) {
-  eval(.parseCheckPlotInputs("ResVsPred"))
+  #----- Validation and formatting of input arguments -----
+  validateIsNotEmpty(data)
+  validateIsOfType(data, "data.frame")
+  dataMapping <- .setDataMapping(dataMapping, ResVsPredDataMapping, data)
   validateIsIncluded(smoother, c("loess", "lm"), nullAllowed = TRUE)
   dataMapping$smoother <- smoother %||% dataMapping$smoother
+
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, ResVsPredPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
+
   mapData <- dataMapping$checkMapData(data)
   mapLabels <- .getAesStringMapping(dataMapping)
 
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
-
-  # Add horizontal lines with offset defined in lines of dataMapping
+  #----- Build layers of molecule plot -----
+  # Each new layer is added on top of previous
+  # Thus, scatter points are added as last layer to prevent them being hidden by lines or errorbars
+  # 1- Diagonal lines
   for (lineIndex in seq_along(dataMapping$lines)) {
-    eval(.parseAddLineLayer("horizontal", dataMapping$lines[[lineIndex]], lineIndex - 1))
-  }
-  if (isEmpty(lineIndex)) {
-    lineIndex <- 0
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = "horizontal",
+      value = dataMapping$lines[[lineIndex]],
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
-  # Add Smoother if defined
+  lineIndex <- ifNotNull(lineIndex, lineIndex, 0)
+
+  # 2- Smoother line
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = lineIndex,
+    plotConfigurationProperty = plotObject$plotConfiguration$lines,
+    propertyNames = c("color", "linetype", "size", "alpha")
+  )
   if (isIncluded(dataMapping$smoother, "loess")) {
     plotObject <- plotObject +
       ggplot2::geom_smooth(
@@ -54,10 +76,10 @@ plotResVsPred <- function(data,
         formula = "y ~ x",
         na.rm = TRUE,
         mapping = ggplot2::aes_string(x = mapLabels$x, y = mapLabels$y),
-        color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-        linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+        color = aestheticValues$color,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        size = aestheticValues$size
       )
   }
   if (isIncluded(dataMapping$smoother, "lm")) {
@@ -69,26 +91,27 @@ plotResVsPred <- function(data,
         se = FALSE,
         formula = "y ~ x",
         na.rm = TRUE,
-        color = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$color, position = lineIndex, aesthetic = "color"),
-        linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$linetype, position = lineIndex, aesthetic = "linetype"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = lineIndex, aesthetic = "alpha"),
-        size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = lineIndex, aesthetic = "size")
+        color = aestheticValues$color,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        size = aestheticValues$size
       )
   }
+  # 2- Scatter points
+  plotObject <- .addScatterLayer(plotObject, data = mapData, mapLabels = mapLabels)
 
-  # If uncertainty is defined, add error bars
-  if (!isOfLength(dataMapping$uncertainty, 0)) {
-    eval(.parseAddUncertaintyLayer())
-  }
-  eval(.parseAddScatterLayer())
-  # Define shapes and colors based on plotConfiguration$points properties
-  eval(.parseUpdateAestheticProperty(AestheticProperties$color, "points"))
-  eval(.parseUpdateAestheticProperty(AestheticProperties$shape, "points"))
-  eval(.parseUpdateAxes())
+  #----- Update properties using ggplot2::scale functions -----
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "points",
+    propertyNames = c("color", "shape"),
+    data = mapData,
+    mapLabels = mapLabels
+  )
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
-
 #' @title plotResVsTime
 #' @description
 #' Producing residuals vs time plots
diff --git a/R/plot-timeprofile.R b/R/plot-timeprofile.R
index 2a6c5773..1b9ca2a1 100644
--- a/R/plot-timeprofile.R
+++ b/R/plot-timeprofile.R
@@ -40,37 +40,47 @@ plotTimeProfile <- function(data = NULL,
                             observedDataMapping = NULL,
                             plotConfiguration = NULL,
                             plotObject = NULL) {
-  validateIsOfType(data, "data.frame", nullAllowed = TRUE)
-  validateIsOfType(observedData, "data.frame", nullAllowed = TRUE)
+  #----- Validation and formatting of input arguments -----
   if (all(isEmpty(data), isEmpty(observedData))) {
-    warning("'data' and 'observedData' are of length 0. Time profile layer was not added.")
-    return(plotObject)
+    stop("At least 'data' or 'observedData' is required.")
   }
+  validateIsOfType(data, "data.frame", nullAllowed = TRUE)
+  validateIsOfType(observedData, "data.frame", nullAllowed = TRUE)
 
   if (!isEmpty(data)) {
-    dataMapping <- dataMapping %||% TimeProfileDataMapping$new(data = data)
+    dataMapping <- .setDataMapping(dataMapping, TimeProfileDataMapping, data)
   }
   if (!isEmpty(observedData)) {
-    observedDataMapping <- observedDataMapping %||% ObservedDataMapping$new(data = data)
+    observedDataMapping <- .setDataMapping(observedDataMapping, ObservedDataMapping, observedData)
   }
+  # If data is empty, plotConfiguration from observedData is used
+  if (isEmpty(data)) {
+    plotConfiguration <- .setPlotConfiguration(
+      plotConfiguration, TimeProfilePlotConfiguration,
+      observedData, metaData, observedDataMapping
+    )
+  }
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, TimeProfilePlotConfiguration,
+    data, metaData, dataMapping
+  )
 
-  plotConfiguration <- plotConfiguration %||% TimeProfilePlotConfiguration$new(data = data, metaData = metaData, dataMapping = dataMapping)
-
-  validateIsOfType(dataMapping, TimeProfileDataMapping, nullAllowed = TRUE)
-  validateIsOfType(observedDataMapping, ObservedDataMapping, nullAllowed = TRUE)
-  validateIsOfType(plotConfiguration, TimeProfilePlotConfiguration)
-
-  # Initialize plot based on plotConfiguration
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
 
-  # Get transformed data from mapping and convert labels into characters usable by aes_string
+  #----- Build layers of molecule plot -----
+  #--- Simulated data ---
   if (!isEmpty(data)) {
     mapData <- dataMapping$checkMapData(data)
     mapLabels <- .getAesStringMapping(dataMapping)
-    # Initialize variables used in legend caption
-    fillValues <- NULL
-    linetypeValues <- NULL
+
+    # Add ribbons for population time profiles
     if (!any(isEmpty(dataMapping$ymin), isEmpty(dataMapping$ymax))) {
+      aestheticValues <- .getAestheticValuesFromConfiguration(
+        n = 1,
+        position = 0,
+        plotConfigurationProperty = plotObject$plotConfiguration$ribbons,
+        propertyNames = c("alpha")
+      )
       plotObject <- plotObject +
         ggplot2::geom_ribbon(
           data = mapData,
@@ -80,11 +90,19 @@ plotTimeProfile <- function(data = NULL,
             ymax = mapLabels$ymax,
             fill = mapLabels$fill
           ),
-          alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$alpha, position = 0, aesthetic = "alpha"),
+          alpha = aestheticValues$alpha,
+          na.rm = TRUE,
           show.legend = TRUE
         )
     }
+    # Add simulated time profile
     if (!isEmpty(dataMapping$y)) {
+      aestheticValues <- .getAestheticValuesFromConfiguration(
+        n = 1,
+        position = 0,
+        plotConfigurationProperty = plotObject$plotConfiguration$lines,
+        propertyNames = c("alpha", "size")
+      )
       plotObject <- plotObject +
         ggplot2::geom_path(
           data = mapData,
@@ -94,215 +112,238 @@ plotTimeProfile <- function(data = NULL,
             color = mapLabels$color,
             linetype = mapLabels$linetype
           ),
-          size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$size, position = 0, aesthetic = "size"),
-          alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$lines$alpha, position = 0, aesthetic = "alpha"),
-          show.legend = TRUE
+          size = aestheticValues$size,
+          alpha = aestheticValues$alpha,
+          na.rm = TRUE,
+          show.legend = TRUE,
         )
     }
-    eval(.parseUpdateAestheticProperty(AestheticProperties$fill, "ribbons"))
-    eval(.parseUpdateAestheticProperty(AestheticProperties$linetype, "lines"))
-    # fillLength defined in .parseUpdateAestheticProperty
-    fillValues <- .getAestheticValues(
-      n = fillLength,
-      selectionKey = plotConfiguration$ribbons$fill,
-      aesthetic = "fill"
+
+    #----- Update properties using ggplot2::scale functions -----
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "ribbons",
+      propertyNames = "fill",
+      data = mapData,
+      mapLabels = mapLabels
+    )
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "lines",
+      propertyNames = "linetype",
+      data = mapData,
+      mapLabels = mapLabels
     )
-    linetypeValues <- .getAestheticValues(
-      n = linetypeLength,
-      selectionKey = plotConfiguration$lines$linetype,
-      aesthetic = "linetype"
+
+    # Get aesthetic information of simulated data for creating a legend caption data.frame
+    simColumnNames <- .getAesPropertyColumnNameFromLabels(
+      mapLabels,
+      c("color", "fill", "linetype")
+    )
+    simAesLengths <- .getAesPropertyLengthFromLabels(
+      data = mapData,
+      simColumnNames,
+      c("color", "fill", "linetype")
     )
+
+    simFillValues <- .getAestheticValuesFromConfiguration(
+      n = simAesLengths$fill,
+      plotConfigurationProperty = plotObject$plotConfiguration$ribbons,
+      propertyNames = "fill"
+    )$fill
+    simLinetypeValues <- .getAestheticValuesFromConfiguration(
+      n = simAesLengths$linetype,
+      plotConfigurationProperty = plotObject$plotConfiguration$lines,
+      propertyNames = "linetype"
+    )$linetype
+    simColorValues <- .getAestheticValuesFromConfiguration(
+      n = simAesLengths$color,
+      plotConfigurationProperty = plotObject$plotConfiguration$lines,
+      propertyNames = "color"
+    )$color
   }
 
-  # If no observed data, also update colors and return plotObect
+  #--- Simulated data only ---
   if (isEmpty(observedData)) {
-    eval(.parseUpdateAestheticProperty(AestheticProperties$color, "lines"))
-    eval(.parseUpdateAxes())
-    # Update and match legend caption to properties
-    # colorLength defined in .parseUpdateAestheticProperty
-    colorValues <- .getAestheticValues(
-      n = colorLength,
-      selectionKey = plotConfiguration$lines$color,
-      aesthetic = "color"
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "lines",
+      propertyNames = "color",
+      data = mapData,
+      mapLabels = mapLabels
     )
-    plotObject$plotConfiguration$legend$caption <- data.frame(
-      name = levels(mapData[, colorVariable]),
-      label = levels(mapData[, colorVariable]),
-      color = colorValues,
-      fill = fillValues %||% NA,
-      linetype = linetypeValues %||% "blank",
-      shape = " ",
-      stringsAsFactors = FALSE
+
+    # Add a legend caption list if legend or properties need to be updated
+    plotObject$plotConfiguration$legend$caption <- list(
+      color = data.frame(
+        names = levels(mapData[, simColumnNames$color]),
+        labels = levels(mapData[, simColumnNames$color]),
+        values = simColorValues,
+        stringsAsFactors = FALSE
+      ),
+      fill = data.frame(
+        names = levels(mapData[, simColumnNames$fill]),
+        labels = levels(mapData[, simColumnNames$fill]),
+        values = simFillValues %||% NA,
+        stringsAsFactors = FALSE
+      ),
+      linetype = data.frame(
+        names = levels(mapData[, simColumnNames$linetype]),
+        labels = levels(mapData[, simColumnNames$linetype]),
+        values = simLinetypeValues %||% "blank",
+        stringsAsFactors = FALSE
+      ),
+      shape = data.frame(
+        names = levels(mapData[, simColumnNames$color]),
+        labels = levels(mapData[, simColumnNames$color]),
+        values = " ",
+        stringsAsFactors = FALSE
+      )
     )
-    eval(.parseUpdateAxes())
+
+    plotObject <- .updateAxes(plotObject)
     return(plotObject)
   }
 
+  #--- Observed data ---
+  # Then, add observed data
   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,
-        mapping = ggplot2::aes_string(
-          x = observedMapLabels$x,
-          ymin = observedMapLabels$ymin,
-          ymax = observedMapLabels$y,
-          color = observedMapLabels$color,
-          group = 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"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, position = 0, aesthetic = "alpha"),
-        show.legend = FALSE
-      ) +
-      ggplot2::geom_linerange(
-        data = mapObservedData,
-        mapping = ggplot2::aes_string(
-          x = observedMapLabels$x,
-          ymin = observedMapLabels$y,
-          ymax = observedMapLabels$ymax,
-          color = observedMapLabels$color,
-          group = 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"),
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, position = 0, aesthetic = "alpha"),
-        show.legend = FALSE
-      ) + 
-      ggplot2::geom_point(
-        data = mapObservedData,
-        mapping = ggplot2::aes_string(
-          x = observedMapLabels$x,
-          y = observedMapLabels$ymin,
-          color = observedMapLabels$color,
-          group = observedMapLabels$color
-        ),
-        size = tlfEnv$defaultErrorbarCapSize,
-        shape = "_",
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, position = 0, aesthetic = "alpha"),
-        show.legend = FALSE
-      ) +
-      ggplot2::geom_point(
-        data = mapObservedData,
-        mapping = ggplot2::aes_string(
-          x = observedMapLabels$x,
-          y = observedMapLabels$ymax,
-          color = observedMapLabels$color,
-          group = observedMapLabels$color
-        ),
-        size = tlfEnv$defaultErrorbarCapSize,
-        shape = "_",
-        alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$errorbars$alpha, position = 0, aesthetic = "alpha"),
-        show.legend = FALSE
-      )
-  }
-  plotObject <- plotObject +
-    ggplot2::geom_point(
+    plotObject <- .addErrorbarLayer(
+      plotObject,
       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"),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$alpha, position = 0, aesthetic = "alpha"),
-      show.legend = TRUE
+      mapLabels = observedMapLabels
     )
+  }
 
-  # 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"
+  plotObject <- .addScatterLayer(
+    plotObject,
+    data = mapObservedData,
+    mapLabels = observedMapLabels
   )
-  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")
+  plotObject <- .updateAesProperties(
+    plotObject,
+    plotConfigurationProperty = "points",
+    propertyNames = "shape",
+    data = mapObservedData,
+    mapLabels = observedMapLabels
+  )
 
-  # No simulated data -> update only observedData
+  obsColumnNames <- .getAesPropertyColumnNameFromLabels(
+    observedMapLabels,
+    c("color", "shape")
+  )
+  obsAesLengths <- .getAesPropertyLengthFromLabels(
+    data = mapObservedData,
+    obsColumnNames,
+    c("color", "shape")
+  )
+
+  obsShapeValues <- .getAestheticValuesFromConfiguration(
+    n = obsAesLengths$shape,
+    plotConfigurationProperty = plotObject$plotConfiguration$points,
+    propertyNames = "shape"
+  )$shape
+  obsColorValues <- .getAestheticValuesFromConfiguration(
+    n = obsAesLengths$color,
+    plotConfigurationProperty = plotObject$plotConfiguration$points,
+    propertyNames = "color"
+  )$color
+
+  #--- Observed data only ---
   if (isEmpty(data)) {
-    colorVariable <- gsub("`", "", observedMapLabels$color)
-    colorLength <- length(levels(mapObservedData[, colorVariable]))
-    colorValues <- .getAestheticValues(
-      n = colorLength,
-      selectionKey = plotConfiguration$points$color,
-      aesthetic = "color"
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "points",
+      propertyNames = "color",
+      data = mapObservedData,
+      mapLabels = observedMapLabels
     )
 
-    suppressMessages(
-      plotObject <- plotObject + ggplot2::scale_color_manual(values = colorValues)
-    )
-    if (isIncluded(colorVariable, "legendLabels")) {
-      plotObject <- plotObject + ggplot2::guides(color = "none")
-    }
-    # Update and match legend caption to properties
-    plotObject$plotConfiguration$legend$caption <- data.frame(
-      name = levels(mapObservedData[, colorVariable]),
-      label = levels(mapObservedData[, colorVariable]),
-      color = colorValues,
-      fill = NA,
-      linetype = "blank",
-      shape = shapeValues %||% NA,
-      stringsAsFactors = FALSE
+    plotObject$plotConfiguration$legend$caption <- list(
+      color = data.frame(
+        names = levels(mapObservedData[, obsColumnNames$color]),
+        labels = levels(mapObservedData[, obsColumnNames$color]),
+        values = obsColorValues,
+        stringsAsFactors = FALSE
+      ),
+      fill = data.frame(
+        names = levels(mapObservedData[, obsColumnNames$color]),
+        labels = levels(mapObservedData[, obsColumnNames$color]),
+        values = NA,
+        stringsAsFactors = FALSE
+      ),
+      linetype = data.frame(
+        names = levels(mapObservedData[, obsColumnNames$color]),
+        labels = levels(mapObservedData[, obsColumnNames$color]),
+        values = "blank",
+        stringsAsFactors = FALSE
+      ),
+      shape = data.frame(
+        names = levels(mapObservedData[, obsColumnNames$shape]),
+        labels = levels(mapObservedData[, obsColumnNames$shape]),
+        values = obsShapeValues %||% NA,
+        stringsAsFactors = FALSE
+      )
     )
 
-    eval(.parseUpdateAxes())
+    plotObject <- .updateAxes(plotObject)
     return(plotObject)
   }
 
-  # Simulated and Observed data -> need to merge the legends
-  colorVariable <- gsub("`", "", mapLabels$color)
-  colorLength <- length(levels(mapData[, colorVariable]))
-  colorObservedVariable <- gsub("`", "", observedMapLabels$color)
+  #--- Simulated and observed data ---
   # 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)
+  colorBreaks <- levels(c(mapData[, simColumnNames$color], mapObservedData[, obsColumnNames$color]))
+  totalColorLength <- length(colorBreaks)
 
   # colorValues are selected colors for simulated (and shared observed data) and then colors for remaining observed data
   # the function ".getAestheticValues" selects these values as defined in the plotConfiguration object
-  colorValues <- c(
+  finalColorValues <- c(
     .getAestheticValues(
-      n = colorLength,
+      n = simAesLengths$color,
       selectionKey = plotConfiguration$lines$color,
       aesthetic = "color"
     ),
-    # By using position = colorLength,
+    # By using position = simAesLengths$color,
     # the function will start selecting the colors that come after the colors selected for lines
     # this aims at preventing a reset of the colors and a need for manual update of the user
     .getAestheticValues(
-      n = totalLength - colorLength,
+      n = totalColorLength - simAesLengths$color,
       selectionKey = plotConfiguration$points$color,
-      position = colorLength,
+      position = simAesLengths$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(fillValues, rep(NA, totalLength - fillLength)),
-    linetype = c(linetypeValues, rep("blank", totalLength - linetypeLength)),
-    shape = c(rep(" ", totalLength - shapeLength), shapeValues),
-    stringsAsFactors = FALSE
+  plotObject$plotConfiguration$legend$caption <- list(
+    color = data.frame(
+      names = colorBreaks,
+      labels = colorBreaks,
+      values = finalColorValues,
+      stringsAsFactors = FALSE
+    ),
+    fill = data.frame(
+      names = levels(mapData[, simColumnNames$fill]),
+      labels = levels(mapData[, simColumnNames$fill]),
+      values = simFillValues %||% NA,
+      stringsAsFactors = FALSE
+    ),
+    linetype = data.frame(
+      names = levels(mapData[, simColumnNames$linetype]),
+      labels = levels(mapData[, simColumnNames$linetype]),
+      values = simLinetypeValues %||% "blank",
+      stringsAsFactors = FALSE
+    ),
+    shape = data.frame(
+      names = levels(mapObservedData[, obsColumnNames$shape]),
+      labels = levels(mapObservedData[, obsColumnNames$shape]),
+      values = obsShapeValues %||% NA,
+      stringsAsFactors = FALSE
+    )
   )
 
   plotObject <- updateTimeProfileLegend(
@@ -310,11 +351,11 @@ plotTimeProfile <- function(data = NULL,
     caption = plotObject$plotConfiguration$legend$caption
   )
 
-  if (isIncluded(colorVariable, "legendLabels") & isIncluded(colorObservedVariable, "legendLabels")) {
+  # remove the legend of aesthetic if default unmapped aesthetic
+  if (isIncluded(simColumnNames$color, "legendLabels") & isIncluded(obsColumnNames$color, "legendLabels")) {
     plotObject <- plotObject + ggplot2::guides(color = "none")
   }
-
-  eval(.parseUpdateAxes())
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
 
@@ -326,39 +367,55 @@ plotTimeProfile <- function(data = NULL,
 #' @return A `ggplot` object
 #' @export
 updateTimeProfileLegend <- function(plotObject, caption) {
-  # Update defined aesthetic properies
-  captionLinetype <- caption[caption$linetype != "blank", ]
-  captionShape <- caption[caption$shape != " ", ]
-  captionFill <- caption[!is.na(caption$fill), ]
+  # Update defined aesthetic properties
+  isDefinedLinetype <- caption$linetype$values != "blank"
+  isDefinedShape <- caption$shape$values != " "
+  isDefinedFill <- !is.na(caption$fill$values)
 
-  if (!isEmpty(captionLinetype)) {
+  if (any(isDefinedLinetype)) {
     suppressMessages(
       plotObject <- plotObject +
-        ggplot2::scale_linetype_manual(breaks = captionLinetype$name, labels = captionLinetype$label, values = captionLinetype$linetype)
+        ggplot2::scale_linetype_manual(
+          breaks = caption$linetype$names[isDefinedLinetype],
+          labels = caption$linetype$labels[isDefinedLinetype],
+          values = caption$linetype$values[isDefinedLinetype]
+        )
     )
   }
-  if (!isEmpty(captionShape)) {
+  if (any(isDefinedShape)) {
     suppressMessages(
       plotObject <- plotObject +
-        ggplot2::scale_shape_manual(breaks = captionShape$name, labels = captionShape$label, values = captionShape$shape)
+        ggplot2::scale_shape_manual(
+          breaks = caption$shape$names[isDefinedShape],
+          labels = caption$shape$labels[isDefinedShape],
+          values = caption$shape$values[isDefinedShape]
+        )
     )
   }
-  if (!isEmpty(captionFill)) {
+  if (any(isDefinedFill)) {
     suppressMessages(
       plotObject <- plotObject +
-        ggplot2::scale_fill_manual(breaks = captionFill$name, labels = captionFill$label, values = captionFill$fill)
+        ggplot2::scale_fill_manual(
+          breaks = caption$fill$names[isDefinedFill],
+          labels = caption$fill$labels[isDefinedFill],
+          values = caption$fill$values[isDefinedFill]
+        )
     )
   }
 
   suppressMessages(
     plotObject <- plotObject +
-      ggplot2::scale_color_manual(breaks = caption$name, labels = caption$label, values = caption$color) +
+      ggplot2::scale_color_manual(
+        breaks = caption$color$names,
+        labels = caption$color$labels,
+        values = caption$color$values
+      ) +
       ggplot2::guides(
         fill = "none", shape = "none", linetype = "none",
         color = ggplot2::guide_legend(
           title = plotObject$plotConfiguration$legend$title$text,
           title.theme = plotObject$plotConfiguration$legend$title$createPlotFont(),
-          override.aes = list(fill = caption$fill, linetype = caption$linetype, shape = caption$shape)
+          override.aes = list(fill = caption$fill$values, linetype = caption$linetype$values)
         )
       )
   )
diff --git a/R/plot-tornado.R b/R/plot-tornado.R
index 9d897565..40882dc1 100644
--- a/R/plot-tornado.R
+++ b/R/plot-tornado.R
@@ -35,11 +35,7 @@ plotTornado <- function(data = NULL,
                         dataMapping = NULL,
                         plotConfiguration = NULL,
                         plotObject = NULL) {
-  validateIsOfType(data, "data.frame", nullAllowed = TRUE)
-  validateIsString(colorPalette, nullAllowed = TRUE)
-  validateIsLogical(sorted, nullAllowed = TRUE)
-  validateIsLogical(bar)
-
+  #----- Validation and formatting of input arguments -----
   if (is.null(data)) {
     validateIsNumeric(x)
     y <- y %||% rep("", length(x))
@@ -53,26 +49,30 @@ plotTornado <- function(data = NULL,
       data = data
     )
   }
+  validateIsNotEmpty(data)
+  validateIsOfType(data, "data.frame")
+  dataMapping <- .setDataMapping(dataMapping, TornadoDataMapping, data)
 
-  dataMapping <- dataMapping %||% TornadoDataMapping$new(data = data)
+  # Update dataMapping if inputs provided by user
+  validateIsLogical(sorted, nullAllowed = TRUE)
   dataMapping$sorted <- sorted %||% dataMapping$sorted
-  plotConfiguration <- plotConfiguration %||% TornadoPlotConfiguration$new(
-    bar = bar,
-    colorPalette = colorPalette,
-    data = data,
-    metaData = metaData,
-    dataMapping = dataMapping
-  )
 
-  validateIsOfType(dataMapping, "TornadoDataMapping")
-  validateIsOfType(plotConfiguration, "TornadoPlotConfiguration")
+  plotConfiguration <- .setPlotConfiguration(
+    plotConfiguration, TornadoPlotConfiguration,
+    data, metaData, dataMapping
+  )
+  # Update plotConfiguration if inputs provided by user
+  validateIsString(colorPalette, nullAllowed = TRUE)
+  validateIsLogical(bar)
+  plotConfiguration$colorPalette <- colorPalette %||% plotConfiguration$colorPalette
+  plotConfiguration$bar <- bar
 
-  plotObject <- plotObject %||% initializePlot(plotConfiguration)
+  plotObject <- .setPlotObject(plotObject, plotConfiguration)
 
-  # Get transformed data from mapping and convert labels into characters usable by aes_string
   mapData <- dataMapping$checkMapData(data)
   mapLabels <- .getAesStringMapping(dataMapping)
 
+  #----- Build layers of molecule plot -----
   # Option sorting the values, which put the wider spread at the top and smaller at the bottom
   # An additional options can be added to change the type of sort for later versions
   # (e.g. increasing/decreasing absolute values or actual values...)
@@ -83,6 +83,13 @@ plotTornado <- function(data = NULL,
   # If tornado is a bar plot, plot configuration option "bar" is TRUE.
   # Otherwise, the plot will use points instead
   if (plotConfiguration$bar) {
+    aestheticValues <- .getAestheticValuesFromConfiguration(
+      n = 1,
+      position = 0,
+      plotConfigurationProperty = plotObject$plotConfiguration$ribbons,
+      propertyNames = c("linetype", "size", "alpha")
+    )
+
     plotObject <- plotObject + ggplot2::geom_col(
       data = mapData,
       mapping = ggplot2::aes_string(
@@ -91,23 +98,30 @@ plotTornado <- function(data = NULL,
         fill = mapLabels$fill,
         color = mapLabels$color
       ),
-      alpha = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$alpha, position = 0, aesthetic = "alpha"),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$size, position = 0, aesthetic = "size"),
-      linetype = .getAestheticValues(n = 1, selectionKey = plotConfiguration$ribbons$linetype, position = 0, aesthetic = "linetype"),
-      position = ggplot2::position_dodge(width = plotConfiguration$dodge)
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size,
+      linetype = aestheticValues$linetype,
+      position = ggplot2::position_dodge(width = plotConfiguration$dodge),
+      na.rm = TRUE
     )
 
-    # Define shapes and colors based on plotConfiguration$points properties
-    fillVariable <- gsub("`", "", mapLabels$fill)
-    colorVariable <- gsub("`", "", mapLabels$color)
-    fillLength <- length(unique(mapData[, fillVariable]))
-    colorLength <- length(unique(mapData[, colorVariable]))
-
-    plotObject <- plotObject +
-      ggplot2::scale_fill_manual(values = .getAestheticValues(n = fillLength, selectionKey = plotConfiguration$ribbons$fill, aesthetic = "fill")) +
-      ggplot2::scale_color_manual(values = .getAestheticValues(n = colorLength, selectionKey = plotConfiguration$ribbons$color, aesthetic = "color"))
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "ribbons",
+      propertyNames = c("fill", "color"),
+      data = mapData,
+      mapLabels = mapLabels
+    )
   }
+
+  # If tornado is a scatter plot, plot configuration option "bar" is FALSE.
   if (!plotConfiguration$bar) {
+    aestheticValues <- .getAestheticValuesFromConfiguration(
+      n = 1,
+      position = 0,
+      plotConfigurationProperty = plotObject$plotConfiguration$points,
+      propertyNames = c("size", "alpha")
+    )
     # For tornado with points, their shape will be taken from the theme properties
     plotObject <- plotObject + ggplot2::geom_point(
       data = mapData,
@@ -117,34 +131,34 @@ plotTornado <- function(data = NULL,
         color = mapLabels$color,
         shape = mapLabels$shape
       ),
-      size = .getAestheticValues(n = 1, selectionKey = plotConfiguration$points$size, position = 0, aesthetic = "size"),
+      size = aestheticValues$size,
+      alpha = aestheticValues$alpha,
       position = ggplot2::position_dodge(width = plotConfiguration$dodge)
     )
 
-    # Define shapes and colors based on plotConfiguration$points properties
-    shapeVariable <- gsub("`", "", mapLabels$shape)
-    colorVariable <- gsub("`", "", mapLabels$color)
-    shapeLength <- length(unique(mapData[, shapeVariable]))
-    colorLength <- length(unique(mapData[, colorVariable]))
-
-    plotObject <- plotObject +
-      ggplot2::scale_shape_manual(values = .getAestheticValues(n = shapeLength, selectionKey = plotConfiguration$points$shape, aesthetic = "shape")) +
-      ggplot2::scale_color_manual(values = .getAestheticValues(n = colorLength, selectionKey = plotConfiguration$points$color, aesthetic = "color"))
+    plotObject <- .updateAesProperties(
+      plotObject,
+      plotConfigurationProperty = "points",
+      propertyNames = c("color", "shape"),
+      data = mapData,
+      mapLabels = mapLabels
+    )
   }
 
-  # Final plot includes a vertical line in 0
-  # And optional color palette otherwise use colors from theme
-  if (!isOfLength(dataMapping$lines, 0)) {
-    plotObject <- plotObject +
-      ggplot2::geom_vline(
-        xintercept = dataMapping$lines,
-        color = .getAestheticValues(n = length(dataMapping$lines), selectionKey = plotConfiguration$lines$color, position = 0, aesthetic = "color"),
-        size = .getAestheticValues(n = length(dataMapping$lines), selectionKey = plotConfiguration$lines$size, position = 0, aesthetic = "size"),
-        linetype = .getAestheticValues(n = length(dataMapping$lines), selectionKey = plotConfiguration$lines$linetype, position = 0, aesthetic = "linetype")
-      )
+  # Add vertical lines
+  for (lineIndex in seq_along(dataMapping$lines)) {
+    plotObject <- .addLineLayer(
+      plotObject,
+      type = "vertical",
+      value = dataMapping$lines[[lineIndex]],
+      # position corresponds to the number of line layers already added
+      position = lineIndex - 1
+    )
   }
 
-  if (!isOfLength(plotConfiguration$colorPalette, 0)) {
+  #----- Update properties using ggplot2::scale functions -----
+  # And optional color palette otherwise use colors from theme
+  if (!isEmpty(plotConfiguration$colorPalette)) {
     try(suppressMessages(
       plotObject <- plotObject +
         ggplot2::scale_fill_brewer(
@@ -153,7 +167,6 @@ plotTornado <- function(data = NULL,
         )
     ))
   }
-  try(suppressMessages(plotObject <- setXAxis(plotObject)))
-  try(suppressMessages(plotObject <- setYAxis(plotObject)))
+  plotObject <- .updateAxes(plotObject)
   return(plotObject)
 }
diff --git a/R/plotconfiguration-legend.R b/R/plotconfiguration-legend.R
index 42d88f51..8c7e7067 100644
--- a/R/plotconfiguration-legend.R
+++ b/R/plotconfiguration-legend.R
@@ -79,7 +79,7 @@ LegendConfiguration <- R6::R6Class(
       if (missing(value)) {
         return(private$.caption)
       }
-      validateIsOfType(value, "data.frame")
+      validateIsOfType(value, c("data.frame", "list"))
       private$.caption <- value
       return(invisible())
     },
diff --git a/R/timeprofile-datamapping.R b/R/timeprofile-datamapping.R
index b5279018..c71f6a30 100644
--- a/R/timeprofile-datamapping.R
+++ b/R/timeprofile-datamapping.R
@@ -32,7 +32,7 @@ TimeProfileDataMapping <- R6::R6Class(
         x = x %||% smartMap$x,
         ymin = ymin %||% smartMap$ymin,
         ymax = ymax %||% smartMap$ymax,
-        color = color, fill = fill, linetype = linetype, group = group
+        color = color, fill = fill, linetype = linetype, group = group, data = data
       )
       # Since TimeProfileDataMapping inherits from RangeDataMapping
       # super$initialize introduce a self$y which is NULL
diff --git a/R/tlf-env.R b/R/tlf-env.R
index 26838e55..38495cbe 100644
--- a/R/tlf-env.R
+++ b/R/tlf-env.R
@@ -189,4 +189,4 @@ setDefaultErrorbarCapSize <- function(size) {
   validateIsNumeric(size)
   tlfEnv$defaultErrorbarCapSize <- size
   return(invisible())
-}
\ No newline at end of file
+}
diff --git a/R/utilities-aesthetics.R b/R/utilities-aesthetics.R
index 1e198c0a..bcc38ef6 100644
--- a/R/utilities-aesthetics.R
+++ b/R/utilities-aesthetics.R
@@ -114,3 +114,101 @@
 .getFirstAestheticValues <- function(n, map) {
   return(.getSameAestheticValues(n, position = 0, map))
 }
+
+
+#' @title .getAesPropertyColumnNameFromLabels
+#' @description Get the column names of the variables mapped to aesthetic properties
+#' @param mapLabels List of mapped label names passed to `ggplot2::aes_string`
+#' @param propertyNames Names of aesthetic property (e.g. `"color"`, `"shape"`...)
+#' @return A list of variable names
+#' @keywords internal
+.getAesPropertyColumnNameFromLabels <- function(mapLabels, propertyNames) {
+  variableNames <- lapply(
+    propertyNames,
+    function(propertyName) {
+      # For aes_string to work with special characters, column names are wrapped with "`"
+      return(gsub("`", "", mapLabels[[propertyName]]))
+    }
+  )
+  names(variableNames) <- propertyNames
+  return(variableNames)
+}
+
+#' @title .getAesPropertyLengthFromLabels
+#' @description Get the names of the variables mapped to aesthetic properties
+#' @param data A data.frame with labels mapped to properties and obtained from a `DataMapping` object
+#' @param columnNames List of mapped column names of `data` obtained
+#' @param propertyNames Names of aesthetic property (e.g. `"color"`, `"shape"`...)
+#' @return A list of variable names
+#' @keywords internal
+.getAesPropertyLengthFromLabels <- function(data, columnNames, propertyNames) {
+  variableLengths <- lapply(
+    propertyNames,
+    function(propertyName) {
+      return(length(unique(data[, columnNames[[propertyName]]])))
+    }
+  )
+  names(variableLengths) <- propertyNames
+  return(variableLengths)
+}
+
+
+#' @title .updateAesProperties
+#' @description Updates the aesthetic properties of `plotObject`
+#' @param plotObject A `ggplot` object
+#' @param plotConfigurationProperty `PlotConfiguration` property name included in .
+#' `"points"`, `"lines"`, `"ribbons"` or `"errorbars"`
+#' @param propertyNames Names of aesthetic property (e.g. `"color"`, `"shape"`...)
+#' @param data A data.frame with labels mapped to properties and obtained from a `DataMapping` object
+#' @param mapLabels List of mapped label names passed to `ggplot2::aes_string`
+#' @return A `ggplot` object
+#' @keywords internal
+#' @import ggplot2
+.updateAesProperties <- function(plotObject, plotConfigurationProperty, propertyNames, data, mapLabels) {
+  propertyVariables <- .getAesPropertyColumnNameFromLabels(mapLabels, propertyNames)
+  propertyLengths <- .getAesPropertyLengthFromLabels(data, propertyVariables, propertyNames)
+  for (propertyName in propertyNames) {
+    # For match.fun to work, ggplot2 namespace is required
+    ggplotScaleFunction <- match.fun(paste0("scale_", propertyName, "_manual"))
+    # Use suppress messages to remove warning of overwriting property scale
+    suppressMessages(
+      plotObject <- plotObject +
+        ggplotScaleFunction(
+          values = .getAestheticValues(
+            n = propertyLengths[[propertyName]],
+            selectionKey = plotObject$plotConfiguration[[plotConfigurationProperty]][[propertyName]],
+            aesthetic = propertyName
+          )
+        )
+    )
+    # remove the legend of aesthetic if default unmapped aesthetic
+    if (isIncluded(propertyVariables[[propertyName]], "legendLabels")) {
+      # Dynamic code is needed here for selecting the correct input parameter of function guides
+      eval(parse(text = paste0(
+        "plotObject <- plotObject + ggplot2::guides(", propertyName, " = 'none')"
+      )))
+    }
+  }
+  return(plotObject)
+}
+
+#' @title .getAestheticValuesFromConfiguration
+#' @description Get list of values for requested aesthetic property
+#' @param n integer defining size of returned aesthetic vector
+#' @param position integer defining the current position in the aesthetic map
+#' @param plotObject A `ggplot` object
+#' @param plotConfigurationProperty `PlotConfiguration` property name included in .
+#' `"points"`, `"lines"`, `"ribbons"` or `"errorbars"`
+#' @param propertyNames Names of aesthetic property (e.g. `"color"`, `"shape"`...)
+#' @return A list of values for requested aesthetic property
+#' @keywords internal
+.getAestheticValuesFromConfiguration <- function(n = 1, position = 0, plotConfigurationProperty, propertyNames) {
+  aestheticValues <- lapply(
+    propertyNames,
+    function(propertyName) {
+      .getAestheticValues(n = n, selectionKey = plotConfigurationProperty[[propertyName]], position = position, aesthetic = propertyName)
+    }
+  )
+  names(aestheticValues) <- propertyNames
+  return(aestheticValues)
+}
diff --git a/R/utilities-molecule-plots.R b/R/utilities-molecule-plots.R
new file mode 100644
index 00000000..a5cf6049
--- /dev/null
+++ b/R/utilities-molecule-plots.R
@@ -0,0 +1,244 @@
+#' @title .addScatterLayer
+#' @description Add scatter points layer of a molecule plot
+#' @param plotObject A `ggplot` object
+#' @param data A data.frame with labels mapped to properties and obtained from a `DataMapping` object
+#' @param mapLabels List of mapped label names passed to `ggplot2::aes_string`
+#' @return A `ggplot` object
+#' @keywords internal
+.addScatterLayer <- function(plotObject, data, mapLabels) {
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = 0,
+    plotConfigurationProperty = plotObject$plotConfiguration$points,
+    propertyNames = c("size", "alpha")
+  )
+
+  plotObject <- plotObject +
+    ggplot2::geom_point(
+      data = data,
+      mapping = ggplot2::aes_string(
+        x = mapLabels$x,
+        y = mapLabels$y,
+        color = mapLabels$color,
+        shape = mapLabels$shape
+      ),
+      size = aestheticValues$size,
+      alpha = aestheticValues$alpha,
+      na.rm = TRUE,
+      show.legend = TRUE
+    )
+  return(plotObject)
+}
+
+#' @title .addErrorbarLayer
+#' @description Add errorbar layer of a molecule plot
+#' @param plotObject A `ggplot` object
+#' @param data A data.frame with labels mapped to properties and obtained from a `DataMapping` object
+#' @param mapLabels List of mapped label names passed to `ggplot2::aes_string`
+#' @return A `ggplot` object
+#' @keywords internal
+.addErrorbarLayer <- function(plotObject, data, mapLabels, direction = "vertical") {
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = 0,
+    plotConfigurationProperty = plotObject$plotConfiguration$errorbars,
+    propertyNames = c("size", "linetype", "alpha")
+  )
+
+  plotObject <- switch(direction,
+    "vertical" = plotObject +
+      ggplot2::geom_linerange(
+        data = data,
+        mapping = ggplot2::aes_string(
+          x = mapLabels$x,
+          # If lower value is negative and plot is log scaled,
+          # Upper bar will still be plotted
+          ymin = mapLabels$ymin,
+          ymax = mapLabels$y,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = aestheticValues$size,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      ggplot2::geom_linerange(
+        data = data,
+        mapping = ggplot2::aes_string(
+          x = mapLabels$x,
+          # If lower value is negative and plot is log scaled,
+          # Upper bar will still be plotted
+          ymin = mapLabels$y,
+          ymax = mapLabels$ymax,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = aestheticValues$size,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      # Add lower cap to error bar
+      # If lower value is negative and plot is log scaled,
+      # Upper bar cap will still be plotted
+      ggplot2::geom_point(
+        data = data,
+        mapping = aes_string(
+          x = mapLabels$x,
+          y = mapLabels$ymin,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = tlfEnv$defaultErrorbarCapSize,
+        shape = "_",
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      # Add upper cap to error bar
+      # If lower value is negative and plot is log scaled,
+      # Upper bar cap will still be plotted
+      ggplot2::geom_point(
+        data = data,
+        mapping = aes_string(
+          x = mapLabels$x,
+          y = mapLabels$ymax,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = tlfEnv$defaultErrorbarCapSize,
+        shape = "_",
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ),
+    "horizontal" = plotObject +
+      ggplot2::geom_linerange(
+        data = data,
+        mapping = ggplot2::aes_string(
+          # If lower value is negative and plot is log scaled,
+          # Upper bar will still be plotted
+          xmin = mapLabels$xmin,
+          xmax = mapLabels$x,
+          y = mapLabels$y,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = aestheticValues$size,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      ggplot2::geom_linerange(
+        data = data,
+        mapping = ggplot2::aes_string(
+          # If lower value is negative and plot is log scaled,
+          # Upper bar will still be plotted
+          xmin = mapLabels$x,
+          xmax = mapLabels$xmax,
+          y = mapLabels$y,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = aestheticValues$size,
+        linetype = aestheticValues$linetype,
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      # Add lower cap to error bar
+      # If lower value is negative and plot is log scaled,
+      # Upper bar cap will still be plotted
+      ggplot2::geom_point(
+        data = data,
+        mapping = aes_string(
+          x = mapLabels$xmin,
+          y = mapLabels$y,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = tlfEnv$defaultErrorbarCapSize,
+        shape = "|",
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      ) +
+      # Add upper cap to error bar
+      # If lower value is negative and plot is log scaled,
+      # Upper bar cap will still be plotted
+      ggplot2::geom_point(
+        data = data,
+        mapping = aes_string(
+          x = mapLabels$xmax,
+          y = mapLabels$y,
+          color = mapLabels$color,
+          group = mapLabels$shape
+        ),
+        size = tlfEnv$defaultErrorbarCapSize,
+        shape = "|",
+        alpha = aestheticValues$alpha,
+        na.rm = TRUE,
+        show.legend = FALSE
+      )
+  )
+  return(plotObject)
+}
+
+#' @title .addLineLayer
+#' @description Add line layer of a molecule plot
+#' @param plotObject A `ggplot` object
+#' @param type one of "horizontal", "vertical" or "diagonal"
+#' Note that for "diagonal", geom_abline is used.
+#' `value` of intercept is taken as is for linear scale but corresponds to the log of `value` for log scale.
+#'  For instance, intercept = c(-1, 0, 1) with log scale actually means that the line will go through c(0.1, 1, 10)
+#'  because c(-1, 0, 1) = log10(c(0.1, 1, 10)).
+#' @param value value of xintercept or yintercept
+#' @param position line position for aesthetic properties
+#' @return A `ggplot` object
+#' @keywords internal
+.addLineLayer <- function(plotObject, type, value, position) {
+  aestheticValues <- .getAestheticValuesFromConfiguration(
+    n = 1,
+    position = position,
+    plotConfigurationProperty = plotObject$plotConfiguration$lines,
+    propertyNames = c("color", "linetype", "size", "alpha")
+  )
+
+  plotObject <- plotObject + switch(type,
+    "horizontal" = ggplot2::geom_hline(
+      yintercept = value,
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
+    ),
+    "vertical" = ggplot2::geom_vline(
+      xintercept = value,
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
+    ),
+    "diagonal" = ggplot2::geom_abline(
+      slope = 1,
+      intercept = value,
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
+    ),
+    "ddiHorizontal" = ggplot2::geom_abline(
+      slope = 0,
+      intercept = value,
+      color = aestheticValues$color,
+      linetype = aestheticValues$linetype,
+      alpha = aestheticValues$alpha,
+      size = aestheticValues$size
+    )
+  )
+  return(plotObject)
+}
diff --git a/man/BoxWhiskerDataMapping.Rd b/man/BoxWhiskerDataMapping.Rd
index 7fb6ba7a..65ac3429 100644
--- a/man/BoxWhiskerDataMapping.Rd
+++ b/man/BoxWhiskerDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{BoxWhiskerDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{BoxWhiskerDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/BoxWhiskerPlotConfiguration.Rd b/man/BoxWhiskerPlotConfiguration.Rd
index 915b7e49..cac5e548 100644
--- a/man/BoxWhiskerPlotConfiguration.Rd
+++ b/man/BoxWhiskerPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{BoxWhiskerPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{BoxWhiskerPlotConfiguration}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/DDIRatioDataMapping.Rd b/man/DDIRatioDataMapping.Rd
index d3745286..148e38b9 100644
--- a/man/DDIRatioDataMapping.Rd
+++ b/man/DDIRatioDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{\link[tlf:PKRatioDataMapping]{tlf::PKRatioDataMapping}} -> \code{DDIRatioDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{tlf::PKRatioDataMapping} -> \code{DDIRatioDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/DDIRatioPlotConfiguration.Rd b/man/DDIRatioPlotConfiguration.Rd
index cf3e1914..3630ca14 100644
--- a/man/DDIRatioPlotConfiguration.Rd
+++ b/man/DDIRatioPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{DDIRatioPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{DDIRatioPlotConfiguration}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/HistogramDataMapping.Rd b/man/HistogramDataMapping.Rd
index 0d7f0759..00d2eec3 100644
--- a/man/HistogramDataMapping.Rd
+++ b/man/HistogramDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{HistogramDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{HistogramDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/HistogramPlotConfiguration.Rd b/man/HistogramPlotConfiguration.Rd
index f2548394..17bb9982 100644
--- a/man/HistogramPlotConfiguration.Rd
+++ b/man/HistogramPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{HistogramPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{HistogramPlotConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/LineElement.Rd b/man/LineElement.Rd
index debf3281..1008fbc5 100644
--- a/man/LineElement.Rd
+++ b/man/LineElement.Rd
@@ -7,7 +7,7 @@
 R6 class defining the properties of background line elements
 }
 \section{Super class}{
-\code{\link[tlf:BackgroundElement]{tlf::BackgroundElement}} -> \code{LineElement}
+\code{tlf::BackgroundElement} -> \code{LineElement}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ObsVsPredDataMapping.Rd b/man/ObsVsPredDataMapping.Rd
index b70466ad..8e73041c 100644
--- a/man/ObsVsPredDataMapping.Rd
+++ b/man/ObsVsPredDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{ObsVsPredDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{ObsVsPredDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
@@ -48,7 +48,6 @@ Other DataMapping classes:
 \item \href{#method-ObsVsPredDataMapping-clone}{\code{ObsVsPredDataMapping$clone()}}
 }
 }
-
 \if{html}{\out{<hr>}}
 \if{html}{\out{<a id="method-ObsVsPredDataMapping-new"></a>}}
 \if{latex}{\out{\hypertarget{method-ObsVsPredDataMapping-new}{}}}
diff --git a/man/ObsVsPredPlotConfiguration.Rd b/man/ObsVsPredPlotConfiguration.Rd
index 9574a573..f86aab69 100644
--- a/man/ObsVsPredPlotConfiguration.Rd
+++ b/man/ObsVsPredPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{ObsVsPredPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{ObsVsPredPlotConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ObservedDataMapping.Rd b/man/ObservedDataMapping.Rd
index e8bb7700..a226854e 100644
--- a/man/ObservedDataMapping.Rd
+++ b/man/ObservedDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{ObservedDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{ObservedDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/PKRatioDataMapping.Rd b/man/PKRatioDataMapping.Rd
index dd55c9bc..a667392f 100644
--- a/man/PKRatioDataMapping.Rd
+++ b/man/PKRatioDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{PKRatioDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{PKRatioDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
@@ -46,7 +46,6 @@ Other DataMapping classes:
 \item \href{#method-PKRatioDataMapping-clone}{\code{PKRatioDataMapping$clone()}}
 }
 }
-
 \if{html}{\out{<hr>}}
 \if{html}{\out{<a id="method-PKRatioDataMapping-new"></a>}}
 \if{latex}{\out{\hypertarget{method-PKRatioDataMapping-new}{}}}
diff --git a/man/PKRatioPlotConfiguration.Rd b/man/PKRatioPlotConfiguration.Rd
index a8a20703..b96233ce 100644
--- a/man/PKRatioPlotConfiguration.Rd
+++ b/man/PKRatioPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{PKRatioPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{PKRatioPlotConfiguration}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/RangeDataMapping.Rd b/man/RangeDataMapping.Rd
index d4322da6..4e688de7 100644
--- a/man/RangeDataMapping.Rd
+++ b/man/RangeDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{RangeDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{RangeDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/ResVsPredDataMapping.Rd b/man/ResVsPredDataMapping.Rd
index 3ef17f26..cc01b404 100644
--- a/man/ResVsPredDataMapping.Rd
+++ b/man/ResVsPredDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{\link[tlf:ObsVsPredDataMapping]{tlf::ObsVsPredDataMapping}} -> \code{ResVsPredDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{tlf::ObsVsPredDataMapping} -> \code{ResVsPredDataMapping}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ResVsPredPlotConfiguration.Rd b/man/ResVsPredPlotConfiguration.Rd
index 1b3f58cf..d04071ec 100644
--- a/man/ResVsPredPlotConfiguration.Rd
+++ b/man/ResVsPredPlotConfiguration.Rd
@@ -7,7 +7,7 @@
 R6 class defining the configuration of a \code{ggplot} object for Res vs Pred/Time plots
 }
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{ResVsPredPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{ResVsPredPlotConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ResVsTimeDataMapping.Rd b/man/ResVsTimeDataMapping.Rd
index fceba4cd..6db8983d 100644
--- a/man/ResVsTimeDataMapping.Rd
+++ b/man/ResVsTimeDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{\link[tlf:ObsVsPredDataMapping]{tlf::ObsVsPredDataMapping}} -> \code{\link[tlf:ResVsPredDataMapping]{tlf::ResVsPredDataMapping}} -> \code{ResVsTimeDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{tlf::ObsVsPredDataMapping} -> \code{tlf::ResVsPredDataMapping} -> \code{ResVsTimeDataMapping}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ResVsTimePlotConfiguration.Rd b/man/ResVsTimePlotConfiguration.Rd
index 537fd32b..ba7b3832 100644
--- a/man/ResVsTimePlotConfiguration.Rd
+++ b/man/ResVsTimePlotConfiguration.Rd
@@ -7,7 +7,7 @@
 R6 class defining the configuration of a \code{ggplot} object for Res vs Pred/Time plots
 }
 \section{Super classes}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{\link[tlf:ResVsPredPlotConfiguration]{tlf::ResVsPredPlotConfiguration}} -> \code{ResVsTimePlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{tlf::ResVsPredPlotConfiguration} -> \code{ResVsTimePlotConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/ThemeAestheticSelections.Rd b/man/ThemeAestheticSelections.Rd
index 55f90af6..071aa959 100644
--- a/man/ThemeAestheticSelections.Rd
+++ b/man/ThemeAestheticSelections.Rd
@@ -7,7 +7,7 @@
 R6 class defining how plot configurations will use aesthetic maps
 }
 \section{Super class}{
-\code{\link[tlf:ThemeAestheticMaps]{tlf::ThemeAestheticMaps}} -> \code{ThemeAestheticSelections}
+\code{tlf::ThemeAestheticMaps} -> \code{ThemeAestheticSelections}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/TimeProfileDataMapping.Rd b/man/TimeProfileDataMapping.Rd
index 14e4973c..05c84477 100644
--- a/man/TimeProfileDataMapping.Rd
+++ b/man/TimeProfileDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{\link[tlf:RangeDataMapping]{tlf::RangeDataMapping}} -> \code{TimeProfileDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{tlf::RangeDataMapping} -> \code{TimeProfileDataMapping}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/TimeProfilePlotConfiguration.Rd b/man/TimeProfilePlotConfiguration.Rd
index 0e509f82..78da917e 100644
--- a/man/TimeProfilePlotConfiguration.Rd
+++ b/man/TimeProfilePlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{TimeProfilePlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{TimeProfilePlotConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/TornadoDataMapping.Rd b/man/TornadoDataMapping.Rd
index 2dce86d1..f6e86421 100644
--- a/man/TornadoDataMapping.Rd
+++ b/man/TornadoDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super classes}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{\link[tlf:XYGDataMapping]{tlf::XYGDataMapping}} -> \code{TornadoDataMapping}
+\code{tlf::XYDataMapping} -> \code{tlf::XYGDataMapping} -> \code{TornadoDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/TornadoPlotConfiguration.Rd b/man/TornadoPlotConfiguration.Rd
index cf765886..f15499b3 100644
--- a/man/TornadoPlotConfiguration.Rd
+++ b/man/TornadoPlotConfiguration.Rd
@@ -18,7 +18,7 @@ Other PlotConfiguration classes:
 }
 \concept{PlotConfiguration classes}
 \section{Super class}{
-\code{\link[tlf:PlotConfiguration]{tlf::PlotConfiguration}} -> \code{TornadoPlotConfiguration}
+\code{tlf::PlotConfiguration} -> \code{TornadoPlotConfiguration}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/XAxisConfiguration.Rd b/man/XAxisConfiguration.Rd
index 208e4b26..9c4e71de 100644
--- a/man/XAxisConfiguration.Rd
+++ b/man/XAxisConfiguration.Rd
@@ -7,7 +7,7 @@
 R6 class defining the configuration of X-axis
 }
 \section{Super class}{
-\code{\link[tlf:AxisConfiguration]{tlf::AxisConfiguration}} -> \code{XAxisConfiguration}
+\code{tlf::AxisConfiguration} -> \code{XAxisConfiguration}
 }
 \section{Methods}{
 \subsection{Public methods}{
diff --git a/man/XYGDataMapping.Rd b/man/XYGDataMapping.Rd
index 50889dd1..ec05e284 100644
--- a/man/XYGDataMapping.Rd
+++ b/man/XYGDataMapping.Rd
@@ -25,7 +25,7 @@ Other DataMapping classes:
 }
 \concept{DataMapping classes}
 \section{Super class}{
-\code{\link[tlf:XYDataMapping]{tlf::XYDataMapping}} -> \code{XYGDataMapping}
+\code{tlf::XYDataMapping} -> \code{XYGDataMapping}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/YAxisConfiguration.Rd b/man/YAxisConfiguration.Rd
index 1fc671d5..bc5e58fc 100644
--- a/man/YAxisConfiguration.Rd
+++ b/man/YAxisConfiguration.Rd
@@ -7,7 +7,7 @@
 R6 class defining the configuration of Y-axis
 }
 \section{Super class}{
-\code{\link[tlf:AxisConfiguration]{tlf::AxisConfiguration}} -> \code{YAxisConfiguration}
+\code{tlf::AxisConfiguration} -> \code{YAxisConfiguration}
 }
 \section{Public fields}{
 \if{html}{\out{<div class="r6-fields">}}
diff --git a/man/dot-addErrorbarLayer.Rd b/man/dot-addErrorbarLayer.Rd
new file mode 100644
index 00000000..7c6b5efb
--- /dev/null
+++ b/man/dot-addErrorbarLayer.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-molecule-plots.R
+\name{.addErrorbarLayer}
+\alias{.addErrorbarLayer}
+\title{.addErrorbarLayer}
+\usage{
+.addErrorbarLayer(plotObject, data, mapLabels, direction = "vertical")
+}
+\arguments{
+\item{plotObject}{A \code{ggplot} object}
+
+\item{data}{A data.frame with labels mapped to properties and obtained from a \code{DataMapping} object}
+
+\item{mapLabels}{List of mapped label names passed to \code{ggplot2::aes_string}}
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Add errorbar layer of a molecule plot
+}
+\keyword{internal}
diff --git a/man/dot-parseAddLineLayer.Rd b/man/dot-addLineLayer.Rd
similarity index 63%
rename from man/dot-parseAddLineLayer.Rd
rename to man/dot-addLineLayer.Rd
index 1130c0a5..955f9d21 100644
--- a/man/dot-parseAddLineLayer.Rd
+++ b/man/dot-addLineLayer.Rd
@@ -1,12 +1,14 @@
 % Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseAddLineLayer}
-\alias{.parseAddLineLayer}
-\title{.parseAddLineLayer}
+% Please edit documentation in R/utilities-molecule-plots.R
+\name{.addLineLayer}
+\alias{.addLineLayer}
+\title{.addLineLayer}
 \usage{
-.parseAddLineLayer(type, value, position)
+.addLineLayer(plotObject, type, value, position)
 }
 \arguments{
+\item{plotObject}{A \code{ggplot} object}
+
 \item{type}{one of "horizontal", "vertical" or "diagonal"
 Note that for "diagonal", geom_abline is used.
 \code{value} of intercept is taken as is for linear scale but corresponds to the log of \code{value} for log scale.
@@ -18,10 +20,9 @@ because c(-1, 0, 1) = log10(c(0.1, 1, 10)).}
 \item{position}{line position for aesthetic properties}
 }
 \value{
-An expression to \code{eval()}
+A \code{ggplot} object
 }
 \description{
-Create an expression that adds scatter plot layer
-TODO: create a vignette explaining how argument \code{lines} in dataMapping is related to this
+Add line layer of a molecule plot
 }
 \keyword{internal}
diff --git a/man/dot-addScatterLayer.Rd b/man/dot-addScatterLayer.Rd
new file mode 100644
index 00000000..f9e47977
--- /dev/null
+++ b/man/dot-addScatterLayer.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-molecule-plots.R
+\name{.addScatterLayer}
+\alias{.addScatterLayer}
+\title{.addScatterLayer}
+\usage{
+.addScatterLayer(plotObject, data, mapLabels)
+}
+\arguments{
+\item{plotObject}{A \code{ggplot} object}
+
+\item{data}{A data.frame with labels mapped to properties and obtained from a \code{DataMapping} object}
+
+\item{mapLabels}{List of mapped label names passed to \code{ggplot2::aes_string}}
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Add scatter points layer of a molecule plot
+}
+\keyword{internal}
diff --git a/man/dot-getAesPropertyColumnNameFromLabels.Rd b/man/dot-getAesPropertyColumnNameFromLabels.Rd
new file mode 100644
index 00000000..3794ed58
--- /dev/null
+++ b/man/dot-getAesPropertyColumnNameFromLabels.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-aesthetics.R
+\name{.getAesPropertyColumnNameFromLabels}
+\alias{.getAesPropertyColumnNameFromLabels}
+\title{.getAesPropertyColumnNameFromLabels}
+\usage{
+.getAesPropertyColumnNameFromLabels(mapLabels, propertyNames)
+}
+\arguments{
+\item{mapLabels}{List of mapped label names passed to \code{ggplot2::aes_string}}
+
+\item{propertyNames}{Names of aesthetic property (e.g. \code{"color"}, \code{"shape"}...)}
+}
+\value{
+A list of variable names
+}
+\description{
+Get the column names of the variables mapped to aesthetic properties
+}
+\keyword{internal}
diff --git a/man/dot-getAesPropertyLengthFromLabels.Rd b/man/dot-getAesPropertyLengthFromLabels.Rd
new file mode 100644
index 00000000..aec818c1
--- /dev/null
+++ b/man/dot-getAesPropertyLengthFromLabels.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-aesthetics.R
+\name{.getAesPropertyLengthFromLabels}
+\alias{.getAesPropertyLengthFromLabels}
+\title{.getAesPropertyLengthFromLabels}
+\usage{
+.getAesPropertyLengthFromLabels(data, columnNames, propertyNames)
+}
+\arguments{
+\item{data}{A data.frame with labels mapped to properties and obtained from a \code{DataMapping} object}
+
+\item{columnNames}{List of mapped column names of \code{data} obtained}
+
+\item{propertyNames}{Names of aesthetic property (e.g. \code{"color"}, \code{"shape"}...)}
+}
+\value{
+A list of variable names
+}
+\description{
+Get the names of the variables mapped to aesthetic properties
+}
+\keyword{internal}
diff --git a/man/dot-getAestheticValuesFromConfiguration.Rd b/man/dot-getAestheticValuesFromConfiguration.Rd
new file mode 100644
index 00000000..50394160
--- /dev/null
+++ b/man/dot-getAestheticValuesFromConfiguration.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-aesthetics.R
+\name{.getAestheticValuesFromConfiguration}
+\alias{.getAestheticValuesFromConfiguration}
+\title{.getAestheticValuesFromConfiguration}
+\usage{
+.getAestheticValuesFromConfiguration(
+  n = 1,
+  position = 0,
+  plotConfigurationProperty,
+  propertyNames
+)
+}
+\arguments{
+\item{n}{integer defining size of returned aesthetic vector}
+
+\item{position}{integer defining the current position in the aesthetic map}
+
+\item{plotConfigurationProperty}{\code{PlotConfiguration} property name included in .
+\code{"points"}, \code{"lines"}, \code{"ribbons"} or \code{"errorbars"}}
+
+\item{propertyNames}{Names of aesthetic property (e.g. \code{"color"}, \code{"shape"}...)}
+
+\item{plotObject}{A \code{ggplot} object}
+}
+\value{
+A list of values for requested aesthetic property
+}
+\description{
+Get list of values for requested aesthetic property
+}
+\keyword{internal}
diff --git a/man/dot-parseAddScatterLayer.Rd b/man/dot-parseAddScatterLayer.Rd
deleted file mode 100644
index 60248c64..00000000
--- a/man/dot-parseAddScatterLayer.Rd
+++ /dev/null
@@ -1,15 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseAddScatterLayer}
-\alias{.parseAddScatterLayer}
-\title{.parseAddScatterLayer}
-\usage{
-.parseAddScatterLayer()
-}
-\value{
-An expression to \code{eval()}
-}
-\description{
-Create an expression that adds scatter plot layer
-}
-\keyword{internal}
diff --git a/man/dot-parseAddUncertaintyLayer.Rd b/man/dot-parseAddUncertaintyLayer.Rd
deleted file mode 100644
index b81cb50a..00000000
--- a/man/dot-parseAddUncertaintyLayer.Rd
+++ /dev/null
@@ -1,16 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseAddUncertaintyLayer}
-\alias{.parseAddUncertaintyLayer}
-\title{.parseAddUncertaintyLayer}
-\usage{
-.parseAddUncertaintyLayer(direction = "vertical")
-}
-\value{
-An expression to \code{eval()}
-}
-\description{
-Create an expression that adds errorbars
-\code{mapLabels} needs to be obtained from \code{DataMapping} objects
-}
-\keyword{internal}
diff --git a/man/dot-parseCheckPlotInputs.Rd b/man/dot-parseCheckPlotInputs.Rd
deleted file mode 100644
index 91c700e5..00000000
--- a/man/dot-parseCheckPlotInputs.Rd
+++ /dev/null
@@ -1,18 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseCheckPlotInputs}
-\alias{.parseCheckPlotInputs}
-\title{.parseCheckPlotInputs}
-\usage{
-.parseCheckPlotInputs(plotType)
-}
-\arguments{
-\item{plotType}{Type of plot (e.g. "PKRatio" for plotPKRatio)}
-}
-\value{
-An expression to \code{eval()}
-}
-\description{
-Create an expression that checks usual plot inputs
-}
-\keyword{internal}
diff --git a/man/dot-parseUpdateAestheticProperty.Rd b/man/dot-parseUpdateAestheticProperty.Rd
deleted file mode 100644
index f417350b..00000000
--- a/man/dot-parseUpdateAestheticProperty.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseUpdateAestheticProperty}
-\alias{.parseUpdateAestheticProperty}
-\title{.parseUpdateAestheticProperty}
-\usage{
-.parseUpdateAestheticProperty(aestheticProperty, plotConfigurationProperty)
-}
-\arguments{
-\item{aestheticProperty}{Name of aesthetic property as defined in \code{AestheticProperties}}
-
-\item{plotConfigurationProperty}{Name of PlotConfiguration property as defined in \code{AestheticProperties}}
-}
-\value{
-An expression to \code{eval()}
-}
-\description{
-Create an expression that updates the aesthetic properties based on
-the information of \code{PlotConfiguration}
-}
-\keyword{internal}
diff --git a/man/dot-parseUpdateAxes.Rd b/man/dot-parseUpdateAxes.Rd
deleted file mode 100644
index effb1e42..00000000
--- a/man/dot-parseUpdateAxes.Rd
+++ /dev/null
@@ -1,15 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/aaa-utilities.R
-\name{.parseUpdateAxes}
-\alias{.parseUpdateAxes}
-\title{.parseUpdateAxes}
-\usage{
-.parseUpdateAxes()
-}
-\value{
-An expression to \code{eval()}
-}
-\description{
-Create an expression that updates the plot axes
-}
-\keyword{internal}
diff --git a/man/dot-setDataMapping.Rd b/man/dot-setDataMapping.Rd
new file mode 100644
index 00000000..b95dd8fc
--- /dev/null
+++ b/man/dot-setDataMapping.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aaa-utilities.R
+\name{.setDataMapping}
+\alias{.setDataMapping}
+\title{.setDataMapping}
+\usage{
+.setDataMapping(dataMapping, DataMappingClass, data = NULL)
+}
+\arguments{
+\item{dataMapping}{A \code{DataMappingClass} object}
+
+\item{DataMappingClass}{Required class for \code{dataMapping}}
+
+\item{data}{A data.frame potentially used for smart mapping}
+}
+\value{
+A \code{DataMapping} object
+}
+\description{
+Set \code{DataMapping} object internally using \code{tlf} default if \code{dataMapping} is not provided
+}
+\keyword{internal}
diff --git a/man/dot-setPlotConfiguration.Rd b/man/dot-setPlotConfiguration.Rd
new file mode 100644
index 00000000..53c1c2f2
--- /dev/null
+++ b/man/dot-setPlotConfiguration.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aaa-utilities.R
+\name{.setPlotConfiguration}
+\alias{.setPlotConfiguration}
+\title{.setPlotConfiguration}
+\usage{
+.setPlotConfiguration(
+  plotConfiguration,
+  PlotConfigurationClass,
+  data = NULL,
+  metaData = NULL,
+  dataMapping = NULL
+)
+}
+\arguments{
+\item{plotConfiguration}{A \code{PlotConfigurationClass} object}
+
+\item{PlotConfigurationClass}{Required class for \code{plotConfiguration}}
+
+\item{data}{A data.frame potentially used for smart plot configuration}
+
+\item{metaData}{A list of meta data potentially used for smart plot configuration}
+
+\item{dataMapping}{A \code{DataMapping} object potentially used for smart plot configuration}
+}
+\value{
+A \code{PlotConfiguration} object
+}
+\description{
+Set \code{PlotConfiguration} object internally using \code{tlf} default if \code{plotConfiguration} is not provided
+}
+\keyword{internal}
diff --git a/man/dot-setPlotObject.Rd b/man/dot-setPlotObject.Rd
new file mode 100644
index 00000000..17307d47
--- /dev/null
+++ b/man/dot-setPlotObject.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aaa-utilities.R
+\name{.setPlotObject}
+\alias{.setPlotObject}
+\title{.setPlotObject}
+\usage{
+.setPlotObject(plotObject, plotConfiguration = NULL)
+}
+\arguments{
+\item{plotObject}{A \code{ggplot} object}
+
+\item{plotConfiguration}{A \code{PlotConfiguration} object}
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Set a \code{ggplot} object associated with its \code{plotConfiguration}
+}
+\keyword{internal}
diff --git a/man/dot-updateAesProperties.Rd b/man/dot-updateAesProperties.Rd
new file mode 100644
index 00000000..55fc67b5
--- /dev/null
+++ b/man/dot-updateAesProperties.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utilities-aesthetics.R
+\name{.updateAesProperties}
+\alias{.updateAesProperties}
+\title{.updateAesProperties}
+\usage{
+.updateAesProperties(
+  plotObject,
+  plotConfigurationProperty,
+  propertyNames,
+  data,
+  mapLabels
+)
+}
+\arguments{
+\item{plotObject}{A \code{ggplot} object}
+
+\item{plotConfigurationProperty}{\code{PlotConfiguration} property name included in .
+\code{"points"}, \code{"lines"}, \code{"ribbons"} or \code{"errorbars"}}
+
+\item{propertyNames}{Names of aesthetic property (e.g. \code{"color"}, \code{"shape"}...)}
+
+\item{data}{A data.frame with labels mapped to properties and obtained from a \code{DataMapping} object}
+
+\item{mapLabels}{List of mapped label names passed to \code{ggplot2::aes_string}}
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Updates the aesthetic properties of \code{plotObject}
+}
+\keyword{internal}
diff --git a/man/dot-updateAxes.Rd b/man/dot-updateAxes.Rd
new file mode 100644
index 00000000..9bc556a4
--- /dev/null
+++ b/man/dot-updateAxes.Rd
@@ -0,0 +1,15 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aaa-utilities.R
+\name{.updateAxes}
+\alias{.updateAxes}
+\title{.updateAxes}
+\usage{
+.updateAxes(plotObject)
+}
+\value{
+A \code{ggplot} object
+}
+\description{
+Updates the plot axes
+}
+\keyword{internal}
diff --git a/vignettes/plot-time-profile.Rmd b/vignettes/plot-time-profile.Rmd
index 0b5d7d87..d67fa078 100644
--- a/vignettes/plot-time-profile.Rmd
+++ b/vignettes/plot-time-profile.Rmd
@@ -537,11 +537,11 @@ To update the plot and have appropriate colors and labels in the legend, it is p
 
 ```{r get legend properties, results='asis'}
 plotLegend <- getLegendCaption(simAndObsTimeProfile)
-knitr::kable(plotLegend)
+knitr::kable(as.data.frame(plotLegend$color))
 ```
 
 ```{r update legend properties}
-plotLegend <- plotLegend[c(1, 3, 2, 4), ]
+plotLegend$color <- plotLegend$color[c(1, 3, 2, 4), ]
 updateTimeProfileLegend(simAndObsTimeProfile, plotLegend)
 ```