diff --git a/petab/v1/math/printer.py b/petab/v1/math/printer.py
index 347c1328..a421989c 100644
--- a/petab/v1/math/printer.py
+++ b/petab/v1/math/printer.py
@@ -37,7 +37,13 @@ def _print_BooleanFalse(self, expr):
     def _print_Pow(self, expr: sp.Pow):
         """Custom printing for the power operator"""
         base, exp = expr.as_base_exp()
-        return f"{self._print(base)} ^ {self._print(exp)}"
+        str_base = self._print(base)
+        str_exp = self._print(exp)
+        if not base.is_Atom:
+            str_base = f"({str_base})"
+        if not exp.is_Atom:
+            str_exp = f"({str_exp})"
+        return f"{str_base} ^ {str_exp}"
 
     def _print_Infinity(self, expr):
         """Custom printing for infinity"""
diff --git a/tests/v1/math/test_math.py b/tests/v1/math/test_math.py
index 940c5340..60bb04b5 100644
--- a/tests/v1/math/test_math.py
+++ b/tests/v1/math/test_math.py
@@ -38,9 +38,11 @@ def test_assumptions():
 
 
 def test_printer():
+    a, b, c, d = sp.symbols("a b c d", real=True)
     assert petab_math_str(None) == ""
     assert petab_math_str(BooleanTrue()) == "true"
     assert petab_math_str(BooleanFalse()) == "false"
+    assert petab_math_str((a + b) ** (c + d)) == "(a + b) ^ (c + d)"
 
 
 def read_cases():