diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py
index 5823233064..4aa34654e2 100644
--- a/python/amici/sbml_import.py
+++ b/python/amici/sbml_import.py
@@ -851,9 +851,15 @@ def computeModelEquationsObjectiveFunction(self, observables=None,
 
         # add user-provided observables or make all species observable
         if(observables):
-            # Replace logX(.) by log(X,.) since symengine cannot parse the former
+            # Replace logX(.) by log(., X) since symengine cannot parse the former
+            # also replace symengine-incompatible sbml log(basis, x)
             for observable in observables:
-                observables[observable]['formula'] = re.sub(r'(^|\W)log(\d+)\(', r'\1log(\2, ', observables[observable]['formula'])
+                observables[observable]['formula'] = re.sub(r'(^|\W)log(\d+)\(', r'\g<1>1/ln(\2)*ln(', observables[observable]['formula'])
+                repl = replaceLogAB(observables[observable]['formula'])
+                if repl != observables[observable]['formula']:
+                    print('Replaced "%s" by "%s", assuming first argument to log() was the basis.' % (observables[observable]['formula'], repl))
+                    observables[observable]['formula'] = repl
+
             self.observables = sp.DenseMatrix([observables[observable]['formula'] for observable in observables])
             observableNames = [observables[observable]['name'] if 'name' in observables[observable].keys()
                                else 'y' + str(index)
@@ -1480,12 +1486,12 @@ class TemplateAmici(Template):
     delimiter = 'TPL_'
 
 
-def assignmentRules2observables(sbml, filter = lambda *_: True):
+def assignmentRules2observables(sbml_model, filter_function=lambda *_: True):
     """Turn assignment rules into observables.
 
     Arguments:
-        sbml: an sbml Model instance
-        filter: callback function taking assignment variable as input and returning True/False to indicate if the respective rule should be turned into an observable
+        sbml_model: an sbml Model instance
+        filter_function: callback function taking assignment variable as input and returning True/False to indicate if the respective rule should be turned into an observable
 
     Returns:
         A dictionary(observableId:{'name':observableNamem,'formula':formulaString})
@@ -1494,14 +1500,14 @@ def assignmentRules2observables(sbml, filter = lambda *_: True):
 
     """
     observables = {}
-    for p in sbml.getListOfParameters():
+    for p in sbml_model.getListOfParameters():
         parameterId = p.getId()
-        if filter(parameterId):
-            observables[parameterId] = {'name': p.getName(), 'formula': sbml.getAssignmentRuleByVariable(parameterId).getFormula()}
+        if filter_function(parameterId):
+            observables[parameterId] = {'name': p.getName(), 'formula': sbml_model.getAssignmentRuleByVariable(parameterId).getFormula()}
 
     for parameterId in observables:
-        sbml.removeRuleByVariable(parameterId)
-        sbml.removeParameter(parameterId)
+        sbml_model.removeRuleByVariable(parameterId)
+        sbml_model.removeParameter(parameterId)
 
     return observables
 
@@ -1545,4 +1551,47 @@ def constantSpeciesToParameters(sbml_model):
             reaction.removeProduct(speciesId)
     
     return transformable
+
+def replaceLogAB(x):
+    """
+    Replace log(a, b) in the given string by ln(b)/ln(a)
+
+    Works for nested parentheses and nested 'log's. This can be used to circumvent
+    the incompatible argument order between symengine (log(x, basis)) and libsbml
+    (log(basis, x)).
+
+    Arguments:
+            x: string to replace
+    Returns:
+            string with replaced 'log's
+    Raises:
+
+    """
+
+    match = re.search(r'(^|\W)log\(', x)
+    if not match:
+        return x
+
+    # index of 'l' of 'log'
+    logStart = match.start() if match.end() - match.start() == 4 else match.start() + 1
+    level = 0 # parenthesis level
+    posComma = -1 # position of comma in log(a,b)
+    for i in range(logStart + 4, len(x)):
+        if x[i] == '(':
+            level += 1
+        elif x[i] == ')':
+            level -= 1
+            if level == -1: break
+        elif x[i] == ',' and level == 0:
+            posComma = i
+
+    if posComma < 0:
+        # was log(a), not log(a,b), so nothing to replace
+        return x
+
+    replacement = '{prefix}ln({a})/ln({basis}){suffix}'.format(prefix = x[:logStart],
+                                                               suffix = x[i+1:],
+                                                               basis = x[logStart+4 : posComma],
+                                                               a = x[posComma+1 : i])
+    return replaceLogAB(replacement)
     
\ No newline at end of file