diff --git a/examples/synthesize_single_neuron_y_direction.py b/examples/synthesize_single_neuron_y_direction.py
new file mode 100644
index 00000000..9c6e6169
--- /dev/null
+++ b/examples/synthesize_single_neuron_y_direction.py
@@ -0,0 +1,90 @@
+# noqa
+
+# Copyright (C) 2021-2024  Blue Brain Project, EPFL
+#
+# SPDX-License-Identifier: Apache-2.0
+
+"""
+Synthesize a single neuron with global direction
+================================================
+
+This example shows how to synthesize a single cell with different y directions
+"""
+
+import json
+from pathlib import Path
+
+import numpy as np
+from morph_tool.transform import rotate
+
+import neurots
+
+
+def run(output_dir, data_dir):
+    """Run the example for generating a single cell."""
+    np.random.seed(42)
+    with open(data_dir / "bio_params.json", encoding="utf-8") as p_file:
+        params = json.load(p_file)
+    # use trunk angle with y_direction awareness
+    params["basal_dendrite"]["orientation"] = {
+        "mode": "pia_constraint",
+        "values": {"form": "step", "params": [1.5, 0.25]},
+    }
+    params["apical_dendrite"]["orientation"] = {
+        "mode": "normal_pia_constraint",
+        "values": {"direction": {"mean": [0.0], "std": [0.0]}},
+    }
+
+    # Initialize a neuron
+    N = neurots.NeuronGrower(
+        input_distributions=data_dir / "bio_distr.json",
+        input_parameters=params,
+    )
+
+    # Grow the neuron
+    neuron = N.grow()
+
+    # Export the synthesized cell
+    neuron.write(output_dir / "generated_cell_orig.asc")
+
+    np.random.seed(42)
+
+    # Initialize a neuron
+    N = neurots.NeuronGrower(
+        input_distributions=data_dir / "bio_distr.json",
+        input_parameters=params,
+        # context={"y_direction": [0.0, 1.0, 0.0]},
+    )
+
+    # Grow the neuron
+    neuron = N.grow()
+
+    # Export the synthesized cell
+    neuron.write(output_dir / "generated_cell_y.asc")
+
+    np.random.seed(42)
+
+    # Initialize a neuron
+    N = neurots.NeuronGrower(
+        input_distributions=data_dir / "bio_distr.json",
+        input_parameters=params,
+        context={"y_direction": [1.0, 0.0, 0.0]},
+    )
+
+    # Grow the neuron
+    neuron = N.grow()
+
+    # Export the synthesized cell
+    neuron.write(output_dir / "generated_cell_x.asc")
+
+    # the rotated neuron should be the same as original one
+
+    rotate(neuron, [[0, -1, 0], [1, 0, 0], [0, 0, 1]])
+    neuron.write(output_dir / "generated_cell_x_rot.asc")
+
+
+if __name__ == "__main__":
+    result_dir = Path("results_single_neuron")
+    result_dir.mkdir(parents=True, exist_ok=True)
+
+    run(result_dir, Path("data"))
diff --git a/neurots/extract_input/from_neurom.py b/neurots/extract_input/from_neurom.py
index 40a29ac8..a318017a 100644
--- a/neurots/extract_input/from_neurom.py
+++ b/neurots/extract_input/from_neurom.py
@@ -9,7 +9,7 @@
 from neurom import stats
 from neurom.features.morphology import trunk_vectors
 
-from neurots.utils import PIA_DIRECTION
+from neurots.utils import Y_DIRECTION
 
 
 def transform_distr(opt_distr):
@@ -130,7 +130,7 @@ def trunk_neurite_3d_angles(pop, neurite_type, bins):
     apical_3d_angles = []
     for morph in pop.morphologies:
         vecs = trunk_vectors(morph, neurite_type=neurite_type)
-        pia_3d_angles += [nm.morphmath.angle_between_vectors(PIA_DIRECTION, vec) for vec in vecs]
+        pia_3d_angles += [nm.morphmath.angle_between_vectors(Y_DIRECTION, vec) for vec in vecs]
         if neurite_type.name != "apical_dendrite":
             apical_ref_vec = trunk_vectors(morph, neurite_type=nm.APICAL_DENDRITE)
             if len(apical_ref_vec) > 0:
diff --git a/neurots/generate/algorithms/abstractgrower.py b/neurots/generate/algorithms/abstractgrower.py
index 2b13cd3b..9cff0ae9 100644
--- a/neurots/generate/algorithms/abstractgrower.py
+++ b/neurots/generate/algorithms/abstractgrower.py
@@ -23,7 +23,7 @@ class AbstractAlgo:
 
     def __init__(self, input_data, params, start_point, context):
         """The TreeGrower Algorithm initialization."""
-        self.context = context
+        self.context = context if context is not None else {}
         self.input_data = copy.deepcopy(input_data)
         self.params = copy.deepcopy(params)
         self.start_point = start_point
diff --git a/neurots/generate/algorithms/basicgrower.py b/neurots/generate/algorithms/basicgrower.py
index 82339faa..e9d752e1 100644
--- a/neurots/generate/algorithms/basicgrower.py
+++ b/neurots/generate/algorithms/basicgrower.py
@@ -52,7 +52,7 @@ def bifurcate(self, current_section):
         Returns:
             tuple[dict, dict]: Two dictionaries containing the two children sections data.
         """
-        dir1, dir2 = self.bif_method()
+        dir1, dir2 = self.bif_method(y_rotation=self.context.get("y_rotation"))
         first_point = np.array(current_section.last_point)
         stop = current_section.stop_criteria
 
diff --git a/neurots/generate/algorithms/tmdgrower.py b/neurots/generate/algorithms/tmdgrower.py
index 9488d0fe..32c7fba1 100644
--- a/neurots/generate/algorithms/tmdgrower.py
+++ b/neurots/generate/algorithms/tmdgrower.py
@@ -161,7 +161,9 @@ def bifurcate(self, current_section):
         self.barcode.remove_bif(current_section.stop_criteria["TMD"].bif_id)
         ang = self.barcode.angles[current_section.stop_criteria["TMD"].bif_id]
 
-        dir1, dir2 = self.bif_method(current_section.history(), angles=ang)
+        dir1, dir2 = self.bif_method(
+            current_section.history(), angles=ang, y_rotation=self.context.get("y_rotation")
+        )
         first_point = np.array(current_section.last_point)
 
         stop1, stop2 = self.get_stop_criteria(current_section)
@@ -249,7 +251,9 @@ def bifurcate(self, current_section):
         first_point = np.array(current_section.last_point)
 
         if current_section.process == "major":
-            dir1, dir2 = bif_methods["directional"](current_section.direction, angles=ang)
+            dir1, dir2 = bif_methods["directional"](
+                current_section.direction, angles=ang, y_rotation=self.context.get("y_rotation")
+            )
 
             if not self._found_last_bif:
                 self.apical_section = current_section.id
@@ -264,7 +268,9 @@ def bifurcate(self, current_section):
                 if not self._found_last_bif:
                     self._found_last_bif = True
         else:
-            dir1, dir2 = self.bif_method(current_section.history(), angles=ang)
+            dir1, dir2 = self.bif_method(
+                current_section.history(), angles=ang, y_rotation=self.context.get("y_rotation")
+            )
             process1 = "secondary"
             process2 = "secondary"
 
diff --git a/neurots/generate/grower.py b/neurots/generate/grower.py
index 71ed4d2d..10df1160 100644
--- a/neurots/generate/grower.py
+++ b/neurots/generate/grower.py
@@ -35,9 +35,11 @@
 from neurots.generate.soma import Soma
 from neurots.generate.soma import SomaGrower
 from neurots.generate.tree import TreeGrower
+from neurots.morphmath import rotation
 from neurots.morphmath import sample
 from neurots.morphmath.utils import normalize_vectors
 from neurots.preprocess import preprocess_inputs
+from neurots.utils import Y_DIRECTION
 from neurots.utils import NeuroTSError
 from neurots.utils import convert_from_legacy_neurite_type
 from neurots.utils import point_to_section_segment
@@ -91,7 +93,7 @@ def __init__(
     ):
         """Constructor of the NeuronGrower class."""
         self.neuron = Morphology()
-        self.context = context
+        self.context = self._process_context(context)
         if rng_or_seed is None or isinstance(
             rng_or_seed, (int, np.integer, SeedSequence, BitGenerator)
         ):
@@ -135,6 +137,20 @@ def __init__(
 
         self._trunk_orientations_class = trunk_orientations_class
 
+    def _process_context(self, context):
+        """Apply some required processing to the context dictionary."""
+        if context is None:
+            return {}
+        if not isinstance(context, dict):
+            return context
+
+        # we ofen need to use the y_direction as a rotation, so we save to it once here
+        if "y_direction" in context:
+            context["y_rotation"] = rotation.rotation_matrix_from_vectors(
+                Y_DIRECTION, context["y_direction"]
+            )
+        return context
+
     def next(self):
         """Call the "next" method of each neurite grower."""
         for grower in list(self.active_neurites):
@@ -363,11 +379,7 @@ def _simple_grow_trunks(self):
                 )
 
     def _3d_angles_grow_trunks(self):
-        """Grow trunk with 3d_angles method via :func:`.orientation.OrientationManager` class.
-
-        Args:
-            input_parameters_with_3d_anglles (dict): input_parameters with fits for 3d angles
-        """
+        """Grow trunk with 3d_angles method via :func:`.orientation.OrientationManager` class."""
         trunk_orientations_manager = self._trunk_orientations_class(
             soma=self.soma_grower.soma,
             parameters=self.input_parameters,
diff --git a/neurots/generate/orientations.py b/neurots/generate/orientations.py
index 76afdeb8..af3d4e11 100644
--- a/neurots/generate/orientations.py
+++ b/neurots/generate/orientations.py
@@ -16,7 +16,7 @@
 from neurots.morphmath import rotation
 from neurots.morphmath import sample
 from neurots.morphmath.utils import normalize_vectors
-from neurots.utils import PIA_DIRECTION
+from neurots.utils import Y_DIRECTION
 from neurots.utils import NeuroTSError
 from neurots.utils import accept_reject
 
@@ -59,7 +59,7 @@ def __init__(self, soma, parameters, distributions, context, rng):
         self._soma = soma
         self._parameters = parameters
         self._distributions = distributions
-        self._context = context
+        self._context = context if context is not None else {}
         self._rng = rng
 
         self._orientations = {}
@@ -222,7 +222,7 @@ def _mode_normal_pia_constraint(self, values_dict, tree_type):
         the second angle to obtain a 3d direction. For multiple apical trees, `mean` and `std`
         should be two lists with lengths equal to number of trees, otherwise it can be a float.
 
-        Pia direction can be overwritten by the parameter 'pia_direction' value.
+        Pia direction can be overwritten by the parameter 'y_direction' value from context.
         """
         n_orientations = sample.n_neurites(self._distributions[tree_type]["num_trees"], self._rng)
         if (
@@ -250,7 +250,7 @@ def _mode_normal_pia_constraint(self, values_dict, tree_type):
 
             phis = self._rng.uniform(0, 2 * np.pi, len(means))
             angles += spherical_angles_to_pia_orientations(
-                phis, thetas, self._parameters.get("pia_direction", None)
+                phis, thetas, self._context.get("y_rotation", None)
             ).tolist()
         return np.array(angles)
 
@@ -260,10 +260,10 @@ def _mode_pia_constraint(self, _, tree_type):
         See :func:`self._sample_trunk_from_3d_angle` for more details on the algorithm.
         """
         n_orientations = sample.n_neurites(self._distributions[tree_type]["num_trees"], self._rng)
-        pia_direction = self._parameters.get("pia_direction", PIA_DIRECTION)
+        y_direction = self._context.get("y_direction", Y_DIRECTION)
         return np.asarray(
             [
-                self._sample_trunk_from_3d_angle(tree_type, pia_direction)
+                self._sample_trunk_from_3d_angle(tree_type, y_direction)
                 for _ in range(n_orientations)
             ]
         )
@@ -302,9 +302,7 @@ def prob(proposal):
             params = self._parameters[tree_type]["orientation"]["values"]["params"]
             p = _prob(val, *params)
 
-            if self._context is not None and self._context.get(
-                "constraints", []
-            ):  # pragma: no cover
+            if self._context.get("constraints", []):  # pragma: no cover
                 for constraint in self._context["constraints"]:
                     if "trunk_prob" in constraint:
                         p *= constraint["trunk_prob"](proposal, self._soma.center)
@@ -477,13 +475,13 @@ def compute_interval_n_tree(soma, n_trees, rng=np.random):
     return phi_intervals, interval_n_trees
 
 
-def spherical_angles_to_pia_orientations(phis, thetas, pia_direction=None):
+def spherical_angles_to_pia_orientations(phis, thetas, y_rotation=None):
     """Compute orientation from spherical angles where thetas are wrt to pia (default=`[0, 1, 0]`).
 
     Args:
         phis (numpy.ndarray): Polar angles.
         thetas (numpy.ndarray): Azimuthal angles.
-        pia_direction (numpy.ndarray): Direction of pia if different from `[0, 1, 0]`.
+        y_rotation (numpy.ndarray): Rotation of y direction if different from `[0, 1, 0]`.
 
     Returns:
         numpy.ndarray: The orientation vectors where each row corresponds to a phi-theta pair.
@@ -492,8 +490,8 @@ def spherical_angles_to_pia_orientations(phis, thetas, pia_direction=None):
     vector = np.column_stack(
         (np.cos(phis) * np.sin(thetas), np.cos(thetas), np.sin(phis) * np.sin(thetas))
     )
-    if pia_direction is not None:
-        vector = vector.dot(rotation.rotation_matrix_from_vectors(PIA_DIRECTION, pia_direction).T)
+    if y_rotation is not None:
+        vector = vector.dot(y_rotation.T)
     return vector
 
 
diff --git a/neurots/generate/tree.py b/neurots/generate/tree.py
index bb2011fb..f85799c4 100644
--- a/neurots/generate/tree.py
+++ b/neurots/generate/tree.py
@@ -165,7 +165,7 @@ def add_section(
         """
         SGrower = section_growers[self.params["metric"]]
         context = copy.deepcopy(self.context)
-        if self.context is not None and "constraints" in self.context:  # pragma: no cover
+        if "constraints" in self.context:  # pragma: no cover
             context["constraints"] = [
                 constraint
                 for constraint in self.context["constraints"]
diff --git a/neurots/morphmath/bifurcation.py b/neurots/morphmath/bifurcation.py
index 75000d81..2e0505ee 100644
--- a/neurots/morphmath/bifurcation.py
+++ b/neurots/morphmath/bifurcation.py
@@ -10,7 +10,7 @@
 from neurots.morphmath.utils import get_random_point
 
 
-def random(random_generator=np.random):
+def random(random_generator=np.random, y_rotation=None):  # pylint: disable=unused-argument
     """Get 3-d coordinates of a new random point.
 
     The distance between the produced point and (0,0,0) is given by the value D.
@@ -18,47 +18,61 @@ def random(random_generator=np.random):
     dir1 = get_random_point(random_generator=random_generator)
     dir2 = get_random_point(random_generator=random_generator)
 
-    return (np.array(dir1), np.array(dir2))
+    return dir1, dir2
 
 
-def symmetric(direction, angles):
+def symmetric(direction, angles, y_rotation=None):
     """Get 3-d coordinates for two new directions at a selected angle."""
-    # phi0 = angles[0] #not used
-    # theta0 = angles[1] #not used
     phi1 = angles[2] / 2.0
     theta1 = angles[3] / 2.0
 
-    dir1 = rt.rotate_vector(direction, [0, 0, 1], phi1)
-    dir1 = rt.rotate_vector(dir1, [1, 0, 0], theta1)
-    dir2 = rt.rotate_vector(direction, [0, 0, 1], -phi1)
-    dir2 = rt.rotate_vector(dir2, [1, 0, 0], -theta1)
+    x_dir = np.array([1, 0, 0])
+    z_dir = np.array([0, 0, 1])
+    if y_rotation is not None:
+        z_dir = y_rotation.dot(z_dir)
+        x_dir = y_rotation.dot(x_dir)
 
-    return (np.array(dir1), np.array(dir2))
+    dir1 = rt.rotate_vector(direction, z_dir, phi1)
+    dir1 = rt.rotate_vector(dir1, x_dir, theta1)
+    dir2 = rt.rotate_vector(direction, z_dir, -phi1)
+    dir2 = rt.rotate_vector(dir2, x_dir, -theta1)
 
+    return dir1, dir2
 
-def bio_oriented(direction, angles):
+
+def bio_oriented(direction, angles, y_rotation=None):
     """Input: init_phi, init_theta, dphi, dtheta."""
     phi0 = angles[0]
     theta0 = angles[1]
     phi1 = angles[2]
     theta1 = angles[3]
 
-    dir1 = rt.rotate_vector(direction, [0, 0, 1], phi0)
-    dir1 = rt.rotate_vector(dir1, [1, 0, 0], theta0)
-    dir2 = rt.rotate_vector(dir1, [0, 0, 1], phi1)
-    dir2 = rt.rotate_vector(dir2, [1, 0, 0], theta1)
+    x_dir = np.array([1, 0, 0])
+    z_dir = np.array([0, 0, 1])
+    if y_rotation is not None:
+        z_dir = y_rotation.dot(z_dir)
+        x_dir = y_rotation.dot(x_dir)
+
+    dir1 = rt.rotate_vector(direction, z_dir, phi0)
+    dir1 = rt.rotate_vector(dir1, x_dir, theta0)
+    dir2 = rt.rotate_vector(dir1, z_dir, phi1)
+    dir2 = rt.rotate_vector(dir2, x_dir, theta1)
 
-    return (np.array(dir1), np.array(dir2))
+    return dir1, dir2
 
 
-def directional(direction, angles):
+def directional(direction, angles, y_rotation=None):
     """Input: init_phi, init_theta, dphi, dtheta."""
-    # phi0 = angles[0] #not used
-    # theta0 = angles[1] #not used
     phi1 = angles[2]
     theta1 = angles[3]
 
-    dir2 = rt.rotate_vector(direction, [0, 0, 1], phi1)
-    dir2 = rt.rotate_vector(dir2, [1, 0, 0], theta1)
+    x_dir = np.array([1, 0, 0])
+    z_dir = np.array([0, 0, 1])
+    if y_rotation is not None:
+        z_dir = y_rotation.dot(z_dir)
+        x_dir = y_rotation.dot(x_dir)
+
+    dir2 = rt.rotate_vector(direction, z_dir, phi1)
+    dir2 = rt.rotate_vector(dir2, x_dir, theta1)
 
-    return (np.array(direction), np.array(dir2))
+    return direction, dir2
diff --git a/neurots/utils.py b/neurots/utils.py
index 45f1b2ca..2d5046b4 100644
--- a/neurots/utils.py
+++ b/neurots/utils.py
@@ -10,7 +10,7 @@
 import numpy as np
 from neurom import COLS
 
-PIA_DIRECTION = [0.0, 1.0, 0.0]
+Y_DIRECTION = [0.0, 1.0, 0.0]
 
 
 class NeuroTSError(Exception):
diff --git a/tests/test_generate_tree.py b/tests/test_generate_tree.py
index 4810ccbe..cb86500c 100644
--- a/tests/test_generate_tree.py
+++ b/tests/test_generate_tree.py
@@ -76,6 +76,21 @@ def test_TreeGrower():
             sections[len(sections) - num - 1],
         )
 
+    grower = NeuronGrower(
+        input_distributions=distributions,
+        input_parameters=params,
+        context={"y_direction": [1, 0, 0]},
+    )
+
+    assert grower.context["y_direction"] == [1, 0, 0]
+    npt.assert_array_equal(
+        grower.context["y_rotation"], np.array([[0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, 0.0, 1.0]])
+    )
+    grower = NeuronGrower(
+        input_distributions=distributions, input_parameters=params, context="OTHER"
+    )
+    assert grower.context == "OTHER"
+
 
 def test_TreeGrower_termination_length():
     np.random.seed(0)
diff --git a/tests/test_morphmath_bifurcation.py b/tests/test_morphmath_bifurcation.py
index 2b74e603..f921046a 100644
--- a/tests/test_morphmath_bifurcation.py
+++ b/tests/test_morphmath_bifurcation.py
@@ -9,6 +9,8 @@
 from numpy.testing import assert_array_almost_equal
 
 from neurots.morphmath import bifurcation as _bf
+from neurots.morphmath import rotation
+from neurots.utils import Y_DIRECTION
 
 
 def test_get_bif_directional():
@@ -31,6 +33,15 @@ def test_get_bif_directional():
             np.array([0.0, 1.0 / np.sqrt(2), 1.0 / np.sqrt(2)]),
         ),
     )
+    rot = rotation.rotation_matrix_from_vectors(Y_DIRECTION, [1, 0, 0]).T
+
+    assert_array_almost_equal(
+        _bf.directional([0.0, 1.0, 0.0], [0.0, 0.0, 0.0, np.pi / 4], y_rotation=rot),
+        (
+            np.array([0.0, 1.0, 0.0]),
+            np.array([0.0, 1.0, 0.0]),
+        ),
+    )
 
 
 def test_get_bif_bio_oriented():
@@ -46,9 +57,20 @@ def test_get_bif_bio_oriented():
         (np.array([0.0, 1.0 / np.sqrt(2), 1.0 / np.sqrt(2)]), np.array([0.0, 0, 1])),
     )
 
+    rot = rotation.rotation_matrix_from_vectors(Y_DIRECTION, [1, 0, 0]).T
+    assert_array_almost_equal(
+        _bf.bio_oriented([0.0, 1.0, 0.0], [0.0, np.pi / 4, 0.0, np.pi / 4], y_rotation=rot),
+        (np.array([0.0, 1.0, 0.0]), np.array([0.0, 1.0, 0.0])),
+    )
+
 
 def test_get_bif_symmetric():
     assert_array_almost_equal(
         _bf.symmetric([0, 0, 1], [1, 1, 1, 1]),
         [[0.0, -0.479426, 0.877583], [0.0, 0.479426, 0.877583]],
     )
+    rot = rotation.rotation_matrix_from_vectors(Y_DIRECTION, [1, 0, 0]).T
+    assert_array_almost_equal(
+        _bf.symmetric([0, 0, 1], [1, 1, 1, 1], y_rotation=rot),
+        [[0.479426, 0.0, 0.877583], [-0.479426, 0.0, 0.877583]],
+    )
diff --git a/tests/test_orientations.py b/tests/test_orientations.py
index 1b415a9a..fb577ecd 100644
--- a/tests/test_orientations.py
+++ b/tests/test_orientations.py
@@ -14,6 +14,8 @@
 
 from neurots.generate import orientations as tested
 from neurots.generate.soma import Soma
+from neurots.morphmath import rotation
+from neurots.utils import Y_DIRECTION
 from neurots.utils import NeuroTSError
 
 
@@ -264,6 +266,29 @@ def test_spherical_angles_to_orientations():
     )
 
 
+def test_spherical_angles_to_pia_orientations():
+    phis = [0.5 * np.pi, np.pi, np.pi]
+
+    thetas = [0.5 * np.pi, np.pi, 0.5 * np.pi]
+
+    expected_orientations = [[0.0, 0.0, 1.0], [0.0, -1.0, 0.0], [-1.0, 0.0, 0.0]]
+
+    npt.assert_allclose(
+        tested.spherical_angles_to_pia_orientations(phis, thetas),
+        expected_orientations,
+        atol=1e-6,
+    )
+
+    expected_orientations = [[0.0, 0.0, 1.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
+    y_rotation = rotation.rotation_matrix_from_vectors(Y_DIRECTION, [1.0, 0.0, 0.0])
+
+    npt.assert_allclose(
+        tested.spherical_angles_to_pia_orientations(phis, thetas, y_rotation=y_rotation),
+        expected_orientations,
+        atol=1e-6,
+    )
+
+
 def test_points_to_orientations():
     origin = np.array([0.0, 0.0, 0.0])
     points = np.array([[2.0, 0.0, 0.0], [0.0, 0.0, 3.0]])
@@ -392,13 +417,13 @@ def test_orientation_manager__mode_normal_pia_constraint():
     expected = np.array([[0.8067661158336028, 0.5513710782140335, 0.21241084824428402]])
     npt.assert_allclose(actual, expected, rtol=1e-5)
 
-    # test other pia direction
-    parameters["pia_direction"] = [1, 0, 0]
+    # test other y direction
+    y_rotation = rotation.rotation_matrix_from_vectors(Y_DIRECTION, [1.0, 0.0, 0.0])
     om = tested.OrientationManager(
         soma=None,
         parameters=parameters,
         distributions=distributions,
-        context=None,
+        context={"y_rotation": y_rotation},
         rng=np.random.default_rng(seed=0),
     )