diff --git a/HISTORY.md b/HISTORY.md
index 381b264a6..2624c5388 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -18,6 +18,15 @@ As long as the above functions are defined correctly, Turing will be able to use
 
 The `Turing.Inference.isgibbscomponent(::MySampler)` interface function still exists, but in this version the default has been changed to `true`, so you should not need to overload this.
 
+## Optimisation interface
+
+The Optim.jl interface has been removed (so you cannot call `Optim.optimize` directly on Turing models).
+You can use the `maximum_likelihood` or `maximum_a_posteriori` functions with an Optim.jl solver instead (via Optimization.jl: see https://docs.sciml.ai/Optimization/stable/optimization_packages/optim/ for documentation of the available solvers).
+
+## Internal changes
+
+The constructors of `OptimLogDensity` have been replaced with a single constructor, `OptimLogDensity(::DynamicPPL.LogDensityFunction)`.
+
 # 0.41.1
 
 The `ModeResult` struct returned by `maximum_a_posteriori` and `maximum_likelihood` can now be wrapped in `InitFromParams()`.
diff --git a/Project.toml b/Project.toml
index 23a8af183..d3042e71a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -41,11 +41,9 @@ StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
 
 [weakdeps]
 DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 
 [extensions]
 TuringDynamicHMCExt = "DynamicHMC"
-TuringOptimExt = ["Optim", "AbstractPPL"]
 
 [compat]
 ADTypes = "1.9"
@@ -72,7 +70,6 @@ LinearAlgebra = "1"
 LogDensityProblems = "2"
 MCMCChains = "5, 6, 7"
 NamedArrays = "0.9, 0.10"
-Optim = "1"
 Optimization = "3, 4, 5"
 OptimizationOptimJL = "0.1, 0.2, 0.3, 0.4"
 OrderedCollections = "1"
@@ -86,7 +83,3 @@ StatsAPI = "1.6"
 StatsBase = "0.32, 0.33, 0.34"
 StatsFuns = "0.8, 0.9, 1"
 julia = "1.10.8"
-
-[extras]
-DynamicHMC = "bbc10e6e-7c05-544b-b16e-64fede858acb"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
diff --git a/ext/TuringOptimExt.jl b/ext/TuringOptimExt.jl
deleted file mode 100644
index ec87467e7..000000000
--- a/ext/TuringOptimExt.jl
+++ /dev/null
@@ -1,198 +0,0 @@
-module TuringOptimExt
-
-using Turing: Turing
-using AbstractPPL: AbstractPPL
-import Turing: DynamicPPL, NamedArrays, Accessors, Optimisation
-using Optim: Optim
-
-####################
-# Optim.jl methods #
-####################
-
-"""
-    Optim.optimize(model::Model, ::MLE, args...; kwargs...)
-
-Compute a maximum likelihood estimate of the `model`.
-
-# Examples
-
-```julia-repl
-@model function f(x)
-    m ~ Normal(0, 1)
-    x ~ Normal(m, 1)
-end
-
-model = f(1.5)
-mle = optimize(model, MLE())
-
-# Use a different optimizer
-mle = optimize(model, MLE(), NelderMead())
-```
-"""
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MLE,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getloglikelihood)
-    init_vals = DynamicPPL.getparams(f.ldf)
-    optimizer = Optim.LBFGS()
-    return _mle_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MLE,
-    init_vals::AbstractArray,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    optimizer = Optim.LBFGS()
-    return _mle_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MLE,
-    optimizer::Optim.AbstractOptimizer,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getloglikelihood)
-    init_vals = DynamicPPL.getparams(f.ldf)
-    return _mle_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MLE,
-    init_vals::AbstractArray,
-    optimizer::Optim.AbstractOptimizer,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    return _mle_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-
-function _mle_optimize(model::DynamicPPL.Model, args...; kwargs...)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getloglikelihood)
-    return _optimize(f, args...; kwargs...)
-end
-
-"""
-    Optim.optimize(model::Model, ::MAP, args...; kwargs...)
-
-Compute a maximum a posterior estimate of the `model`.
-
-# Examples
-
-```julia-repl
-@model function f(x)
-    m ~ Normal(0, 1)
-    x ~ Normal(m, 1)
-end
-
-model = f(1.5)
-map_est = optimize(model, MAP())
-
-# Use a different optimizer
-map_est = optimize(model, MAP(), NelderMead())
-```
-"""
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MAP,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
-    init_vals = DynamicPPL.getparams(f.ldf)
-    optimizer = Optim.LBFGS()
-    return _map_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MAP,
-    init_vals::AbstractArray,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    optimizer = Optim.LBFGS()
-    return _map_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MAP,
-    optimizer::Optim.AbstractOptimizer,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
-    init_vals = DynamicPPL.getparams(f.ldf)
-    return _map_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-function Optim.optimize(
-    model::DynamicPPL.Model,
-    ::Optimisation.MAP,
-    init_vals::AbstractArray,
-    optimizer::Optim.AbstractOptimizer,
-    options::Optim.Options=Optim.Options();
-    kwargs...,
-)
-    return _map_optimize(model, init_vals, optimizer, options; kwargs...)
-end
-
-function _map_optimize(model::DynamicPPL.Model, args...; kwargs...)
-    f = Optimisation.OptimLogDensity(model, DynamicPPL.getlogjoint)
-    return _optimize(f, args...; kwargs...)
-end
-
-"""
-    _optimize(f::OptimLogDensity, optimizer=Optim.LBFGS(), args...; kwargs...)
-
-Estimate a mode, i.e., compute a MLE or MAP estimate.
-"""
-function _optimize(
-    f::Optimisation.OptimLogDensity,
-    init_vals::AbstractArray=DynamicPPL.getparams(f.ldf),
-    optimizer::Optim.AbstractOptimizer=Optim.LBFGS(),
-    options::Optim.Options=Optim.Options(),
-    args...;
-    kwargs...,
-)
-    # Convert the initial values, since it is assumed that users provide them
-    # in the constrained space.
-    # TODO(penelopeysm): As with in src/optimisation/Optimisation.jl, unclear
-    # whether initialisation is really necessary at all
-    vi = DynamicPPL.unflatten(f.ldf.varinfo, init_vals)
-    vi = DynamicPPL.link(vi, f.ldf.model)
-    f = Optimisation.OptimLogDensity(
-        f.ldf.model, f.ldf.getlogdensity, vi; adtype=f.ldf.adtype
-    )
-    init_vals = DynamicPPL.getparams(f.ldf)
-
-    # Optimize!
-    M = Optim.optimize(Optim.only_fg!(f), init_vals, optimizer, options, args...; kwargs...)
-
-    # Warn the user if the optimization did not converge.
-    if !Optim.converged(M)
-        @warn """
-            Optimization did not converge! You may need to correct your model or adjust the
-            Optim parameters.
-        """
-    end
-
-    # Get the optimum in unconstrained space. `getparams` does the invlinking.
-    vi = f.ldf.varinfo
-    vi_optimum = DynamicPPL.unflatten(vi, M.minimizer)
-    logdensity_optimum = Optimisation.OptimLogDensity(
-        f.ldf.model, f.ldf.getlogdensity, vi_optimum; adtype=f.ldf.adtype
-    )
-    vals_dict = Turing.Inference.getparams(f.ldf.model, vi_optimum)
-    iters = map(AbstractPPL.varname_and_value_leaves, keys(vals_dict), values(vals_dict))
-    vns_vals_iter = mapreduce(collect, vcat, iters)
-    varnames = map(Symbol ∘ first, vns_vals_iter)
-    vals = map(last, vns_vals_iter)
-    vmat = NamedArrays.NamedArray(vals, varnames)
-    return Optimisation.ModeResult(vmat, M, -M.minimum, logdensity_optimum, vals_dict)
-end
-
-end # module
diff --git a/src/optimisation/Optimisation.jl b/src/optimisation/Optimisation.jl
index a183619e0..8c956500b 100644
--- a/src/optimisation/Optimisation.jl
+++ b/src/optimisation/Optimisation.jl
@@ -19,137 +19,65 @@ using StatsAPI: StatsAPI
 using Statistics: Statistics
 using LinearAlgebra: LinearAlgebra
 
-export maximum_a_posteriori, maximum_likelihood
-# The MAP and MLE exports are only needed for the Optim.jl interface.
-export MAP, MLE
+export maximum_a_posteriori, maximum_likelihood, MAP, MLE
 
 """
     ModeEstimator
 
 An abstract type to mark whether mode estimation is to be done with maximum a posteriori
-(MAP) or maximum likelihood estimation (MLE). This is only needed for the Optim.jl interface.
+(MAP) or maximum likelihood estimation (MLE).
 """
 abstract type ModeEstimator end
 
 """
     MLE <: ModeEstimator
 
-Concrete type for maximum likelihood estimation. Only used for the Optim.jl interface.
+Concrete type for maximum likelihood estimation.
 """
 struct MLE <: ModeEstimator end
 
 """
     MAP <: ModeEstimator
 
-Concrete type for maximum a posteriori estimation. Only used for the Optim.jl interface.
+Concrete type for maximum a posteriori estimation.
 """
 struct MAP <: ModeEstimator end
 
 """
-    OptimLogDensity{
-        M<:DynamicPPL.Model,
-        F<:Function,
-        V<:DynamicPPL.AbstractVarInfo,
-        AD<:ADTypes.AbstractADType,
-    }
-
-A struct that wraps a single LogDensityFunction. Can be invoked either using
-
-```julia
-OptimLogDensity(model, varinfo; adtype=adtype)
-```
-
-or
-
-```julia
-OptimLogDensity(model; adtype=adtype)
-```
-
-If not specified, `adtype` defaults to `AutoForwardDiff()`.
+    OptimLogDensity{L<:DynamicPPL.LogDensityFunction}
 
-An OptimLogDensity does not, in itself, obey the LogDensityProblems interface.
-Thus, if you want to calculate the log density of its contents at the point
-`z`, you should manually call
+A struct that represents a log-density function, which can be used with Optimization.jl.
+This is a thin wrapper around `DynamicPPL.LogDensityFunction`: the main difference is that
+the log-density is negated (because Optimization.jl performs minimisation, and we usually
+want to maximise the log-density).
 
-```julia
-LogDensityProblems.logdensity(f.ldf, z)
-```
+An `OptimLogDensity` does not, in itself, obey the LogDensityProblems.jl interface. Thus, if
+you want to calculate the log density of its contents at the point `z`, you should manually
+call `LogDensityProblems.logdensity(f.ldf, z)`, instead of `LogDensityProblems.logdensity(f,
+z)`.
 
-However, it is a callable object which returns the *negative* log density of
-the underlying LogDensityFunction at the point `z`. This is done to satisfy
-the Optim.jl interface.
-
-```julia
-optim_ld = OptimLogDensity(model, varinfo)
-optim_ld(z)  # returns -logp
-```
+However, because Optimization.jl requires the objective function to be callable, you can
+also call `f(z)` directly to get the negative log density at `z`.
 """
 struct OptimLogDensity{L<:DynamicPPL.LogDensityFunction}
     ldf::L
-
-    function OptimLogDensity(
-        model::DynamicPPL.Model,
-        getlogdensity::Function,
-        vi::DynamicPPL.AbstractVarInfo;
-        adtype::ADTypes.AbstractADType=Turing.DEFAULT_ADTYPE,
-    )
-        ldf = DynamicPPL.LogDensityFunction(model, getlogdensity, vi; adtype=adtype)
-        return new{typeof(ldf)}(ldf)
-    end
-    function OptimLogDensity(
-        model::DynamicPPL.Model,
-        getlogdensity::Function;
-        adtype::ADTypes.AbstractADType=Turing.DEFAULT_ADTYPE,
-    )
-        # No varinfo
-        return OptimLogDensity(
-            model,
-            getlogdensity,
-            DynamicPPL.ldf_default_varinfo(model, getlogdensity);
-            adtype=adtype,
-        )
-    end
 end
 
 """
     (f::OptimLogDensity)(z)
     (f::OptimLogDensity)(z, _)
 
-Evaluate the negative log joint or log likelihood at the array `z`. Which one is evaluated
-depends on the context of `f`.
+Evaluate the negative log probability density at the array `z`. Which kind of probability
+density is evaluated depends on the `getlogdensity` function used to construct the
+underlying `LogDensityFunction` (e.g., `DynamicPPL.getlogjoint` for MAP estimation, or
+`DynamicPPL.getloglikelihood` for MLE).
 
-Any second argument is ignored. The two-argument method only exists to match interface the
+Any second argument is ignored. The two-argument method only exists to match the interface
 required by Optimization.jl.
 """
 (f::OptimLogDensity)(z::AbstractVector) = -LogDensityProblems.logdensity(f.ldf, z)
 (f::OptimLogDensity)(z, _) = f(z)
 
-# NOTE: The format of this function is dictated by Optim. The first argument sets whether to
-# compute the function value, the second whether to compute the gradient (and stores the
-# gradient). The last one is the actual argument of the objective function.
-function (f::OptimLogDensity)(F, G, z)
-    if G !== nothing
-        # Calculate log joint and its gradient.
-        logp, ∇logp = LogDensityProblems.logdensity_and_gradient(f.ldf, z)
-
-        # Save the negative gradient to the pre-allocated array.
-        copyto!(G, -∇logp)
-
-        # If F is something, the negative log joint is requested as well.
-        # We have already computed it as a by-product above and hence return it directly.
-        if F !== nothing
-            return -logp
-        end
-    end
-
-    # Only negative log joint requested but no gradient.
-    if F !== nothing
-        return -LogDensityProblems.logdensity(f.ldf, z)
-    end
-
-    return nothing
-end
-
 """
     ModeResult{
         V<:NamedArrays.NamedArray,
@@ -299,7 +227,9 @@ function StatsBase.informationmatrix(
     if linked
         new_vi = DynamicPPL.invlink!!(old_ldf.varinfo, old_ldf.model)
         new_f = OptimLogDensity(
-            old_ldf.model, old_ldf.getlogdensity, new_vi; adtype=old_ldf.adtype
+            DynamicPPL.LogDensityFunction(
+                old_ldf.model, old_ldf.getlogdensity, new_vi; adtype=old_ldf.adtype
+            ),
         )
         m = Accessors.@set m.f = new_f
     end
@@ -314,7 +244,12 @@ function StatsBase.informationmatrix(
         invlinked_ldf = m.f.ldf
         new_vi = DynamicPPL.link!!(invlinked_ldf.varinfo, invlinked_ldf.model)
         new_f = OptimLogDensity(
-            invlinked_ldf.model, old_ldf.getlogdensity, new_vi; adtype=invlinked_ldf.adtype
+            DynamicPPL.LogDensityFunction(
+                invlinked_ldf.model,
+                old_ldf.getlogdensity,
+                new_vi;
+                adtype=invlinked_ldf.adtype,
+            ),
         )
         m = Accessors.@set m.f = new_f
     end
@@ -568,7 +503,8 @@ function estimate_mode(
         vi = DynamicPPL.link(vi, model)
     end
 
-    log_density = OptimLogDensity(model, getlogdensity, vi)
+    ldf = DynamicPPL.LogDensityFunction(model, getlogdensity, vi; adtype=adtype)
+    log_density = OptimLogDensity(ldf)
 
     prob = Optimization.OptimizationProblem(log_density, adtype, constraints)
     solution = Optimization.solve(prob, solver; kwargs...)
diff --git a/test/ext/OptimInterface.jl b/test/ext/OptimInterface.jl
deleted file mode 100644
index 721e255f3..000000000
--- a/test/ext/OptimInterface.jl
+++ /dev/null
@@ -1,194 +0,0 @@
-module OptimInterfaceTests
-
-using ..Models: gdemo_default
-using Distributions.FillArrays: Zeros
-using AbstractPPL: AbstractPPL
-using DynamicPPL: DynamicPPL
-using LinearAlgebra: I
-using Optim: Optim
-using Random: Random
-using StatsBase: StatsBase
-using StatsBase: coef, coefnames, coeftable, informationmatrix, stderror, vcov
-using Test: @test, @testset
-using Turing
-
-@testset "TuringOptimExt" begin
-    @testset "MLE" begin
-        Random.seed!(222)
-        true_value = [0.0625, 1.75]
-
-        m1 = Optim.optimize(gdemo_default, MLE())
-        m2 = Optim.optimize(gdemo_default, MLE(), Optim.NelderMead())
-        m3 = Optim.optimize(gdemo_default, MLE(), true_value, Optim.LBFGS())
-        m4 = Optim.optimize(gdemo_default, MLE(), true_value)
-
-        @test all(isapprox.(m1.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m2.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m3.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m4.values.array - true_value, 0.0, atol=0.01))
-    end
-
-    @testset "MAP" begin
-        Random.seed!(222)
-        true_value = [49 / 54, 7 / 6]
-
-        m1 = Optim.optimize(gdemo_default, MAP())
-        m2 = Optim.optimize(gdemo_default, MAP(), Optim.NelderMead())
-        m3 = Optim.optimize(gdemo_default, MAP(), true_value, Optim.LBFGS())
-        m4 = Optim.optimize(gdemo_default, MAP(), true_value)
-
-        @test all(isapprox.(m1.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m2.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m3.values.array - true_value, 0.0, atol=0.01))
-        @test all(isapprox.(m4.values.array - true_value, 0.0, atol=0.01))
-    end
-
-    @testset "StatsBase integration" begin
-        Random.seed!(54321)
-        mle_est = Optim.optimize(gdemo_default, MLE())
-        # Calculated based on the two data points in gdemo_default, [1.5, 2.0]
-        true_values = [0.0625, 1.75]
-
-        @test coefnames(mle_est) == [:s, :m]
-
-        diffs = coef(mle_est).array - [0.0625031; 1.75001]
-        @test all(isapprox.(diffs, 0.0, atol=0.1))
-
-        infomat = [2/(2 * true_values[1]^2) 0.0; 0.0 2/true_values[1]]
-        @test all(isapprox.(infomat - informationmatrix(mle_est), 0.0, atol=0.01))
-
-        vcovmat = [2 * true_values[1]^2/2 0.0; 0.0 true_values[1]/2]
-        @test all(isapprox.(vcovmat - vcov(mle_est), 0.0, atol=0.01))
-
-        ctable = coeftable(mle_est)
-        @test ctable isa StatsBase.CoefTable
-
-        s = stderror(mle_est).array
-        @test all(isapprox.(s - [0.06250415643292194, 0.17677963626053916], 0.0, atol=0.01))
-
-        @test coefnames(mle_est) == Distributions.params(mle_est)
-        @test vcov(mle_est) == inv(informationmatrix(mle_est))
-
-        @test isapprox(loglikelihood(mle_est), -0.0652883561466624, atol=0.01)
-    end
-
-    @testset "Linear regression test" begin
-        @model function regtest(x, y)
-            beta ~ MvNormal(Zeros(2), I)
-            mu = x * beta
-            return y ~ MvNormal(mu, I)
-        end
-
-        Random.seed!(987)
-        true_beta = [1.0, -2.2]
-        x = rand(40, 2)
-        y = x * true_beta
-
-        model = regtest(x, y)
-        mle = Optim.optimize(model, MLE())
-
-        vcmat = inv(x'x)
-        vcmat_mle = vcov(mle).array
-
-        @test isapprox(mle.values.array, true_beta)
-        @test isapprox(vcmat, vcmat_mle)
-    end
-
-    @testset "Dot tilde test" begin
-        @model function dot_gdemo(x)
-            s ~ InverseGamma(2, 3)
-            m ~ Normal(0, sqrt(s))
-
-            return (.~)(x, Normal(m, sqrt(s)))
-        end
-
-        model_dot = dot_gdemo([1.5, 2.0])
-
-        mle1 = Optim.optimize(gdemo_default, MLE())
-        mle2 = Optim.optimize(model_dot, MLE())
-
-        map1 = Optim.optimize(gdemo_default, MAP())
-        map2 = Optim.optimize(model_dot, MAP())
-
-        @test isapprox(mle1.values.array, mle2.values.array)
-        @test isapprox(map1.values.array, map2.values.array)
-    end
-
-    # FIXME: Some models don't work for ReverseDiff.
-    # TODO: Check if above statement is still correct
-    @testset "MAP for $(model.f)" for model in DynamicPPL.TestUtils.DEMO_MODELS
-        result_true = DynamicPPL.TestUtils.posterior_optima(model)
-
-        optimizers = [Optim.LBFGS(), Optim.NelderMead()]
-        @testset "$(nameof(typeof(optimizer)))" for optimizer in optimizers
-            result = Optim.optimize(model, MAP(), optimizer)
-            vals = result.values
-
-            for vn in DynamicPPL.TestUtils.varnames(model)
-                for vn_leaf in AbstractPPL.varname_leaves(vn, get(result_true, vn))
-                    @test get(result_true, vn_leaf) ≈ vals[Symbol(vn_leaf)] atol = 0.05
-                end
-            end
-        end
-    end
-
-    # Some of the models have one variance parameter per observation, and so
-    # the MLE should have the variances set to 0. Since we're working in
-    # transformed space, this corresponds to `-Inf`, which is of course not achievable.
-    # In particular, it can result in "early termination" of the optimization process
-    # because we hit NaNs, etc. To avoid this, we set the `g_tol` and the `f_tol` to
-    # something larger than the default.
-    allowed_incorrect_mle = [
-        DynamicPPL.TestUtils.demo_dot_assume_observe,
-        DynamicPPL.TestUtils.demo_assume_index_observe,
-        DynamicPPL.TestUtils.demo_assume_multivariate_observe,
-        DynamicPPL.TestUtils.demo_assume_multivariate_observe_literal,
-        DynamicPPL.TestUtils.demo_dot_assume_observe_submodel,
-        DynamicPPL.TestUtils.demo_dot_assume_observe_matrix_index,
-        DynamicPPL.TestUtils.demo_assume_submodel_observe_index_literal,
-        DynamicPPL.TestUtils.demo_dot_assume_observe_index,
-        DynamicPPL.TestUtils.demo_dot_assume_observe_index_literal,
-        DynamicPPL.TestUtils.demo_assume_matrix_observe_matrix_index,
-    ]
-    @testset "MLE for $(model.f)" for model in DynamicPPL.TestUtils.DEMO_MODELS
-        result_true = DynamicPPL.TestUtils.likelihood_optima(model)
-
-        # `NelderMead` seems to struggle with convergence here, so we exclude it.
-        @testset "$(nameof(typeof(optimizer)))" for optimizer in [Optim.LBFGS()]
-            options = Optim.Options(; g_tol=1e-3, f_tol=1e-3)
-            result = Optim.optimize(model, MLE(), optimizer, options)
-            vals = result.values
-
-            for vn in DynamicPPL.TestUtils.varnames(model)
-                for vn_leaf in AbstractPPL.varname_leaves(vn, get(result_true, vn))
-                    if model.f in allowed_incorrect_mle
-                        @test isfinite(get(result_true, vn_leaf))
-                    else
-                        @test get(result_true, vn_leaf) ≈ vals[Symbol(vn_leaf)] atol = 0.05
-                    end
-                end
-            end
-        end
-    end
-
-    # Issue: https://discourse.julialang.org/t/turing-mixture-models-with-dirichlet-weightings/112910
-    @testset "with different linked dimensionality" begin
-        @model demo_dirichlet() = x ~ Dirichlet(2 * ones(3))
-        model = demo_dirichlet()
-        result = Optim.optimize(model, MAP())
-        @test result.values ≈ mode(Dirichlet(2 * ones(3))) atol = 0.2
-    end
-
-    @testset "with :=" begin
-        @model function demo_track()
-            x ~ Normal()
-            return y := 100 + x
-        end
-        model = demo_track()
-        result = Optim.optimize(model, MAP())
-        @test result.values[:x] ≈ 0 atol = 1e-1
-        @test result.values[:y] ≈ 100 atol = 1e-1
-    end
-end
-
-end
diff --git a/test/optimisation/Optimisation.jl b/test/optimisation/Optimisation.jl
index 23c8a7d8e..6ef5fe8f1 100644
--- a/test/optimisation/Optimisation.jl
+++ b/test/optimisation/Optimisation.jl
@@ -44,8 +44,12 @@ using Turing
             m2 = model2() | (x=x,)
             # Doesn't matter if we use getlogjoint or getlogjoint_internal since the
             # VarInfo isn't linked.
-            ld1 = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
-            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint_internal)
+            ld1 = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogjoint)
+            )
+            ld2 = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m2, DynamicPPL.getlogjoint_internal)
+            )
             @test ld1(w) == ld2(w)
         end
 
@@ -53,26 +57,36 @@ using Turing
             vn = @varname(inner)
             m1 = prefix(model1(x), vn)
             m2 = prefix((model2() | (x=x,)), vn)
-            ld1 = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
-            ld2 = Turing.Optimisation.OptimLogDensity(m2, DynamicPPL.getlogjoint_internal)
+            ld1 = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogjoint)
+            )
+            ld2 = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m2, DynamicPPL.getlogjoint_internal)
+            )
             @test ld1(w) == ld2(w)
         end
 
         @testset "Joint, prior, and likelihood" begin
             m1 = model1(x)
             a = [0.3]
-            ld_joint = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogjoint)
-            ld_prior = Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogprior)
+            ld_joint = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogjoint)
+            )
+            ld_prior = Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogprior)
+            )
             ld_likelihood = Turing.Optimisation.OptimLogDensity(
-                m1, DynamicPPL.getloglikelihood
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getloglikelihood)
             )
             @test ld_joint(a) == ld_prior(a) + ld_likelihood(a)
 
             # test that the prior accumulator is calculating the right thing
-            @test Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogprior)([0.3]) ≈
-                -Distributions.logpdf(Uniform(0, 2), 0.3)
-            @test Turing.Optimisation.OptimLogDensity(m1, DynamicPPL.getlogprior)([-0.3]) ≈
-                -Distributions.logpdf(Uniform(0, 2), -0.3)
+            @test Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogprior)
+            )([0.3]) ≈ -Distributions.logpdf(Uniform(0, 2), 0.3)
+            @test Turing.Optimisation.OptimLogDensity(
+                DynamicPPL.LogDensityFunction(m1, DynamicPPL.getlogprior)
+            )([-0.3]) ≈ -Distributions.logpdf(Uniform(0, 2), -0.3)
         end
     end
 
@@ -646,7 +660,9 @@ using Turing
             return nothing
         end
         m = saddle_model()
-        optim_ld = Turing.Optimisation.OptimLogDensity(m, DynamicPPL.getloglikelihood)
+        optim_ld = Turing.Optimisation.OptimLogDensity(
+            DynamicPPL.LogDensityFunction(m, DynamicPPL.getloglikelihood)
+        )
         vals = Turing.Optimisation.NamedArrays.NamedArray([0.0, 0.0])
         m = Turing.Optimisation.ModeResult(
             vals,
diff --git a/test/runtests.jl b/test/runtests.jl
index 81b4bdde2..4ca2e38fe 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -68,7 +68,6 @@ end
 
         @testset "mode estimation" verbose = true begin
             @timeit_include("optimisation/Optimisation.jl")
-            @timeit_include("ext/OptimInterface.jl")
         end
     end