diff --git a/NAMESPACE b/NAMESPACE
index ee3ba55f..b6656d68 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(mcstate_model)
+export(mcstate_model_properties)
 export(mcstate_rng)
 export(mcstate_rng_distributed_pointer)
 export(mcstate_rng_distributed_state)
diff --git a/R/model.R b/R/model.R
index d6a71c35..bf26480e 100644
--- a/R/model.R
+++ b/R/model.R
@@ -1,63 +1,226 @@
-##' Create a basic `mcstate` model.  Currently nothing here is
-##' validated, and it's likely that users will never actually use this
-##' directly.  Contains data and methods that define a basic model
-##' object, so that we can implement samplers against.  Not all models
-##' will support everything here, and we'll add additional
-##' fields/traits over time to advertise what a model can do.  For
-##' example, models will need to advertise that they are capable of
-##' being differentiated, or that they are stochastic in order to be
-##' used with different methods.
+##' Describe properties of a model.  Use of this function is optional,
+##' but you can pass the return value of this as the `properties`
+##' argument of mcstate_model to enforce that your model does actually have these
+##' properties.
 ##'
-##' @title Create basic model
+##' @title Describe model properties
 ##'
-##' @param parameters Names of the parameters.  Every parameter is
-##'   named, and for now every parameter is a scalar.  We might relax
-##'   this later to support an `odin`-style structured parameter list,
-##'   but that might just generate a suitable vector of parameter
-##'   names perhaps?  In any case, once we start doing inference it's
-##'   naturally in the R^n, and here n is defined as the length of
-##'   this vector of names.
+##' @param has_gradient Logical, indicating if the model has a
+##'   `gradient` method.  Use `NULL` (the default) to detect this from
+##'   the model.
 ##'
-##' @param direct_sample A function to sample directly from the
-##'   parameter space, given an [mcstate_rng] object to sample from.
-##'   In the case where a model returns a posterior (e.g., in Bayesian
-##'   inference), this is assumed to be sampling from the prior.
-##'   We'll use this for generating initial conditions for MCMC where
-##'   those are not given, and possibly other uses.
-##'
-##' @param density Compute the model density for a vector of parameter
-##'   values; this is the posterior probability in the case of
-##'   Bayesian inference, but it could be anything really.  Models can
-##'   return `-Inf` if things are impossible, and we'll try and cope
-##'   gracefully with that wherever possible.
-##'
-##' @param gradient Compute the gradient of `density` with respect to
-##'   the parameter vector; takes a parameter vector and returns a
-##'   vector the same length.  For efficiency, the model may want to
-##'   be stateful so that gradients can be efficiently calculated
-##'   after a density calculation, or density after gradient, where
-##'   these are called with the same parameters.
-##'
-##' @param domain Information on the parameter domain.  This is a two
+##' @param has_direct_sample Logical, indicating if the model has a
+##'   `direct_sample` method.  Use `NULL` (the default) to detect this
+##'   from the model.
+##'
+##' @return A list of class `mcstate_model_properties` which should
+##'   not be modified.
+##'
+##' @export
+mcstate_model_properties <- function(has_gradient = NULL,
+                                     has_direct_sample = NULL) {
+  ret <- list(has_gradient = has_gradient,
+              has_direct_sample = has_direct_sample)
+  class(ret) <- "mcstate_model_properties"
+  ret
+}
+
+##' Create a basic `mcstate` model.  This takes a user-supplied object
+##' that minimally can compute a probability density (via a `density`
+##' function) and information about parameters; with this we can
+##' sample from the model using `MCMC` using [mcstate_sample].  We
+##' don't imagine that many users will call this function directly,
+##' but that this will be glue used by packages.
+##'
+##' The `model` argument can be a list or environment (something
+##' indexable by `$`) and have elements:
+##'
+##' * `density`: A function that will compute some probability
+##'   density.  It must take an argument representing a parameter
+##'   vector (a numeric vector) and return a single value.  This is
+##'   the posterior probability density in Bayesian inference, but it
+##'   could be anything really.  Models can return `-Inf` if things
+##'   are impossible, and we'll try and cope gracefully with that
+##'   wherever possible.
+##'
+##' * `parameters`: A character vector of parameter names.  This
+##'   vector is the source of truth for the length of the parameter
+##'   vector.
+##'
+##' * `domain`: Information on the parameter domain.  This is a two
 ##'   column matrix with `length(parameters)` rows representing each
 ##'   parameter.  The parameter minimum and maximum bounds are given
 ##'   as the first and second column.  Infinite values (`-Inf` or
 ##'   `Inf`) should be used where the parameter has infinite domain up
 ##'   or down.  Currently used to translate from a bounded to
 ##'   unbounded space for HMC, but we might also use this for
-##'   reflecting proposals in MCMC too.
+##'   reflecting proposals in MCMC too.  If not present we assume that
+##'   the model is valid everywhere (i.e., that all parameters are
+##'   valid from `-Inf` to `Inf`.
 ##'
-##' @return An object of class `mcstate_model`, which can be used with
-##'   a sampler.
+##' * `direct_sample`: A function to sample directly from the
+##'   parameter space, given an [mcstate_rng] object to sample from.
+##'   In the case where a model returns a posterior (e.g., in Bayesian
+##'   inference), this is assumed to be sampling from the prior.
+##'   We'll use this for generating initial conditions for MCMC where
+##'   those are not given, and possibly other uses.  If not given then
+##'   when using [mcstate_sample()] the user will have to provide a
+##'   vector of initial states.
+##'
+##' * `gradient`: A function to compute the gradient of `density` with
+##'   respect to the parameter vector; takes a parameter vector and
+##'   returns a vector the same length.  For efficiency, the model may
+##'   want to be stateful so that gradients can be efficiently
+##'   calculated after a density calculation, or density after
+##'   gradient, where these are called with the same parameters.  This
+##'   function is optional (and may not be well defined or possible to
+##'   define).
+##'
+##' @title Create basic model
+##'
+##' @param model A list or environment with elements as described in
+##'   Details.
+##'
+##' @param properties Optionally, a [mcstate_model_properties] object,
+##'   used to enforce or clarify properties of the model.
+##'
+##' @return An object of class `mcstate_model`.  This will have elements:
+##'
+##' * `model`: The model as provided
+##' * `parameters`: The parameter name vector
+##' * `domain`: The parameter domain matrix
+##' * `direct_sample`: The `direct_sample` function, if provided by the model
+##' * `gradient`: The `gradient` function, if provided by the model
+##' * `properties`: A list of properties of the model;
+##'   see [mcstate_model_properties()].  Currently this contains:
+##'     * `has_gradient`: the model can compute its gradient
+##'     * `has_direct_sample`: the model can sample from parameters space
 ##'
 ##' @export
-mcstate_model <- function(parameters, direct_sample, density, gradient,
-                          domain) {
-  ret <- list(parameters = parameters,
-              direct_sample = direct_sample,
+mcstate_model <- function(model, properties = NULL) {
+  call <- environment() # for nicer stack traces
+  parameters <- validate_model_parameters(model, call)
+  domain <- validate_model_domain(model, call)
+  density <- validate_model_density(model, call)
+
+  properties <- validate_model_properties(properties, call)
+  gradient <- validate_model_gradient(model, properties, call)
+  direct_sample <- validate_model_direct_sample(model, properties, call)
+
+  ## Update properties based on what we found:
+  properties$has_gradient <- !is.null(gradient)
+  properties$has_direct_sample <- !is.null(direct_sample)
+
+  ret <- list(model = model,
+              parameters = parameters,
+              domain = domain,
               density = density,
               gradient = gradient,
-              domain = domain)
+              direct_sample = direct_sample,
+              properties = properties)
   class(ret) <- "mcstate_model"
   ret
 }
+
+
+
+validate_model_properties <- function(properties, call = NULL) {
+  if (is.null(properties)) {
+    return(mcstate_model_properties())
+  }
+
+  if (!inherits(properties, "mcstate_model_properties")) {
+    cli::cli_abort(
+      "Expected 'properties' to be a 'mcstate_model_properties' object",
+      arg = "properties", call = call)
+  }
+
+  properties
+}
+
+
+validate_model_parameters <- function(model, call = NULL) {
+  if (!is.character(model$parameters)) {
+    cli::cli_abort("Expected 'model$parameters' to be a character vector",
+                   arg = "model", call = call)
+  }
+  model$parameters
+}
+
+
+validate_model_domain <- function(model, call = NULL) {
+  n_pars <- length(model$parameters)
+  if (is.null(model$domain)) {
+    domain <- cbind(rep(-Inf, n_pars), rep(Inf, n_pars))
+  } else {
+    domain <- model$domain
+    if (!is.matrix(domain)) {
+      cli::cli_abort("Expected 'model$domain' to be a matrix if non-NULL")
+    }
+    if (nrow(domain) != n_pars) {
+      cli::cli_abort(paste(
+        "Expected 'model$domain' to have {n_pars} row{?s},",
+        "but it had {nrow(domain)}"))
+    }
+    if (ncol(domain) != 2) {
+      cli::cli_abort(paste(
+        "Expected 'model$domain' to have 2 columns,",
+        "but it had {ncol(domain)}"))
+    }
+  }
+  domain
+}
+
+
+validate_model_density <- function(model, call = NULL) {
+  if (!is.function(model$density)) {
+    cli::cli_abort("Expected 'model$density' to be a function",
+                   arg = "model", call = call)
+  }
+  model$density
+}
+
+
+validate_model_optional_method <- function(model, properties,
+                                           method_name, property_name,
+                                           call) {
+  if (isFALSE(properties[[property_name]])) {
+    return(NULL)
+  }
+  value <- model[[method_name]]
+  if (isTRUE(properties[[property_name]]) && !is.function(value)) {
+    cli::cli_abort(
+      paste("Did not find a function '{method_name}' within your model, but",
+            "your properties say that it should do"),
+      arg = "model", call = call)
+  }
+  if (!is.null(value) && !is.function(value)) {
+    cli::cli_abort(
+      "Expected 'model${method_name}' to be a function if non-NULL",
+      arg = "model", call = call)
+  }
+  value
+}
+
+
+validate_model_gradient <- function(model, properties, call) {
+  validate_model_optional_method(
+    model, properties, "gradient", "has_gradient", call)
+}
+
+
+validate_model_direct_sample <- function(model, properties, call) {
+  validate_model_optional_method(
+    model, properties, "direct_sample", "has_direct_sample", call)
+}
+
+
+require_direct_sample <- function(model, message, ...) {
+  if (!model$properties$has_direct_sample) {
+    cli::cli_abort(
+      c(message,
+        i = paste("This model does not provide 'direct_sample()', so we",
+                  "cannot directly sample from its parameter space")),
+      ...)
+  }
+}
diff --git a/R/sample.R b/R/sample.R
index 5b736a3f..983807c9 100644
--- a/R/sample.R
+++ b/R/sample.R
@@ -73,6 +73,9 @@ initial_parameters <- function(initial, model, rng, call = NULL) {
   n_pars <- length(model$parameters)
   n_chains <- length(rng)
   if (is.null(initial)) {
+    require_direct_sample(model,
+                          "'initial' must be provided with this model",
+                          arg = "initial", call = environment())
     ## Really this would just be from the prior; we can't directly
     ## sample from the posterior!
     initial <- lapply(rng, function(r) model$direct_sample(r))
diff --git a/man/mcstate_model.Rd b/man/mcstate_model.Rd
index f5271354..b4832af5 100644
--- a/man/mcstate_model.Rd
+++ b/man/mcstate_model.Rd
@@ -4,58 +4,78 @@
 \alias{mcstate_model}
 \title{Create basic model}
 \usage{
-mcstate_model(parameters, direct_sample, density, gradient, domain)
+mcstate_model(model, properties = NULL)
 }
 \arguments{
-\item{parameters}{Names of the parameters.  Every parameter is
-named, and for now every parameter is a scalar.  We might relax
-this later to support an \code{odin}-style structured parameter list,
-but that might just generate a suitable vector of parameter
-names perhaps?  In any case, once we start doing inference it's
-naturally in the R^n, and here n is defined as the length of
-this vector of names.}
+\item{model}{A list or environment with elements as described in
+Details.}
 
-\item{direct_sample}{A function to sample directly from the
-parameter space, given an \link{mcstate_rng} object to sample from.
-In the case where a model returns a posterior (e.g., in Bayesian
-inference), this is assumed to be sampling from the prior.
-We'll use this for generating initial conditions for MCMC where
-those are not given, and possibly other uses.}
-
-\item{density}{Compute the model density for a vector of parameter
-values; this is the posterior probability in the case of
-Bayesian inference, but it could be anything really.  Models can
-return \code{-Inf} if things are impossible, and we'll try and cope
-gracefully with that wherever possible.}
-
-\item{gradient}{Compute the gradient of \code{density} with respect to
-the parameter vector; takes a parameter vector and returns a
-vector the same length.  For efficiency, the model may want to
-be stateful so that gradients can be efficiently calculated
-after a density calculation, or density after gradient, where
-these are called with the same parameters.}
-
-\item{domain}{Information on the parameter domain.  This is a two
+\item{properties}{Optionally, a \link{mcstate_model_properties} object,
+used to enforce or clarify properties of the model.}
+}
+\value{
+An object of class \code{mcstate_model}.  This will have elements:
+\itemize{
+\item \code{model}: The model as provided
+\item \code{parameters}: The parameter name vector
+\item \code{domain}: The parameter domain matrix
+\item \code{direct_sample}: The \code{direct_sample} function, if provided by the model
+\item \code{gradient}: The \code{gradient} function, if provided by the model
+\item \code{properties}: A list of properties of the model;
+see \code{\link[=mcstate_model_properties]{mcstate_model_properties()}}.  Currently this contains:
+\itemize{
+\item \code{has_gradient}: the model can compute its gradient
+\item \code{has_direct_sample}: the model can sample from parameters space
+}
+}
+}
+\description{
+Create a basic \code{mcstate} model.  This takes a user-supplied object
+that minimally can compute a probability density (via a \code{density}
+function) and information about parameters; with this we can
+sample from the model using \code{MCMC} using \link{mcstate_sample}.  We
+don't imagine that many users will call this function directly,
+but that this will be glue used by packages.
+}
+\details{
+The \code{model} argument can be a list or environment (something
+indexable by \code{$}) and have elements:
+\itemize{
+\item \code{density}: A function that will compute some probability
+density.  It must take an argument representing a parameter
+vector (a numeric vector) and return a single value.  This is
+the posterior probability density in Bayesian inference, but it
+could be anything really.  Models can return \code{-Inf} if things
+are impossible, and we'll try and cope gracefully with that
+wherever possible.
+\item \code{parameters}: A character vector of parameter names.  This
+vector is the source of truth for the length of the parameter
+vector.
+\item \code{domain}: Information on the parameter domain.  This is a two
 column matrix with \code{length(parameters)} rows representing each
 parameter.  The parameter minimum and maximum bounds are given
 as the first and second column.  Infinite values (\code{-Inf} or
 \code{Inf}) should be used where the parameter has infinite domain up
 or down.  Currently used to translate from a bounded to
 unbounded space for HMC, but we might also use this for
-reflecting proposals in MCMC too.}
-}
-\value{
-An object of class \code{mcstate_model}, which can be used with
-a sampler.
+reflecting proposals in MCMC too.  If not present we assume that
+the model is valid everywhere (i.e., that all parameters are
+valid from \code{-Inf} to \code{Inf}.
+\item \code{direct_sample}: A function to sample directly from the
+parameter space, given an \link{mcstate_rng} object to sample from.
+In the case where a model returns a posterior (e.g., in Bayesian
+inference), this is assumed to be sampling from the prior.
+We'll use this for generating initial conditions for MCMC where
+those are not given, and possibly other uses.  If not given then
+when using \code{\link[=mcstate_sample]{mcstate_sample()}} the user will have to provide a
+vector of initial states.
+\item \code{gradient}: A function to compute the gradient of \code{density} with
+respect to the parameter vector; takes a parameter vector and
+returns a vector the same length.  For efficiency, the model may
+want to be stateful so that gradients can be efficiently
+calculated after a density calculation, or density after
+gradient, where these are called with the same parameters.  This
+function is optional (and may not be well defined or possible to
+define).
 }
-\description{
-Create a basic \code{mcstate} model.  Currently nothing here is
-validated, and it's likely that users will never actually use this
-directly.  Contains data and methods that define a basic model
-object, so that we can implement samplers against.  Not all models
-will support everything here, and we'll add additional
-fields/traits over time to advertise what a model can do.  For
-example, models will need to advertise that they are capable of
-being differentiated, or that they are stochastic in order to be
-used with different methods.
 }
diff --git a/man/mcstate_model_properties.Rd b/man/mcstate_model_properties.Rd
new file mode 100644
index 00000000..716168c3
--- /dev/null
+++ b/man/mcstate_model_properties.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model.R
+\name{mcstate_model_properties}
+\alias{mcstate_model_properties}
+\title{Describe model properties}
+\usage{
+mcstate_model_properties(has_gradient = NULL, has_direct_sample = NULL)
+}
+\arguments{
+\item{has_gradient}{Logical, indicating if the model has a
+\code{gradient} method.  Use \code{NULL} (the default) to detect this from
+the model.}
+
+\item{has_direct_sample}{Logical, indicating if the model has a
+\code{direct_sample} method.  Use \code{NULL} (the default) to detect this
+from the model.}
+}
+\value{
+A list of class \code{mcstate_model_properties} which should
+not be modified.
+}
+\description{
+Describe properties of a model.  Use of this function is optional,
+but you can pass the return value of this as the \code{properties}
+argument of mcstate_model to enforce that your model does actually have these
+properties.
+}
diff --git a/tests/testthat/helper-mcstate2.R b/tests/testthat/helper-mcstate2.R
index 36fb58b9..893de0db 100644
--- a/tests/testthat/helper-mcstate2.R
+++ b/tests/testthat/helper-mcstate2.R
@@ -4,12 +4,12 @@ ex_simple_gamma1 <- function(shape = 1, rate = 1) {
   e$rate <- rate
   with(
     e,
-    mcstate_model(
+    mcstate_model(list(
       parameters = "gamma",
       direct_sample = function(rng) {
         rng$gamma(1, shape = shape, scale = 1 / rate)
       },
       density = function(x) dgamma(x, shape = shape, rate = rate, log = TRUE),
       gradient = function(x) (shape - 1) / x - rate,
-      domain = rbind(c(0, Inf))))
+      domain = rbind(c(0, Inf)))))
 }
diff --git a/tests/testthat/test-model.R b/tests/testthat/test-model.R
index 28970a56..743e3622 100644
--- a/tests/testthat/test-model.R
+++ b/tests/testthat/test-model.R
@@ -1,13 +1,119 @@
-test_that("can create a basic model", {
-  m <- local({
-    shape <- 1
-    rate <- 1
-    mcstate_model(
-      parameters = "gamma",
-      direct_sample = function() rgamma(1, shape = shape, rate = rate),
-      density = function(x) dgamma(x, shape = shape, rate = rate, log = TRUE),
-      gradient = function(x) (shape - 1) / x - shape,
-      domain = rbind(c(0, Inf)))
-  })
+test_that("can create a minimal model", {
+  m <- mcstate_model(list(density = function(x) dnorm(x, log = TRUE),
+                          parameters = "a"))
   expect_s3_class(m, "mcstate_model")
+  expect_equal(m$properties,
+               mcstate_model_properties(has_gradient = FALSE,
+                                        has_direct_sample = FALSE))
+  expect_equal(m$domain, cbind(-Inf, Inf))
+  expect_equal(m$parameters, "a")
+  expect_equal(m$density(0), dnorm(0, log = TRUE))
+})
+
+
+test_that("can create a more interesting model", {
+  m <- ex_simple_gamma1()
+  expect_equal(m$properties,
+               mcstate_model_properties(has_gradient = TRUE,
+                                        has_direct_sample = TRUE))
+  expect_equal(m$domain, cbind(0, Inf))
+  expect_equal(m$parameters, "gamma")
+  expect_equal(m$density(1), dgamma(1, 1, 1, log = TRUE))
+})
+
+
+test_that("Require parameters to be given if model does not provide them", {
+  expect_error(
+    mcstate_model(list(density = function(x) dnorm(x, log = TRUE))),
+    "Expected 'model$parameters' to be a character vector",
+    fixed = TRUE)
+})
+
+
+test_that("require density is a function", {
+  expect_error(mcstate_model(list(parameters = "a")),
+               "Expected 'model$density' to be a function",
+               fixed = TRUE)
+})
+
+
+test_that("require gradient is a function if given", {
+  expect_error(
+    mcstate_model(list(density = identity, gradient = TRUE, parameters = "a")),
+    "Expected 'model$gradient' to be a function if non-NULL",
+    fixed = TRUE)
+})
+
+
+test_that("require direct sample is a function if given", {
+  expect_error(
+    mcstate_model(list(density = identity,
+                       direct_sample = TRUE,
+                       parameters = "a")),
+    "Expected 'model$direct_sample' to be a function if non-NULL",
+    fixed = TRUE)
+})
+
+
+test_that("validate domain", {
+  expect_error(
+    mcstate_model(list(density = identity, parameters = "a", domain = list())),
+    "Expected 'model$domain' to be a matrix if non-NULL",
+    fixed = TRUE)
+  expect_error(
+    mcstate_model(list(density = identity,
+                       parameters = "a",
+                       domain = matrix(1:4, 2, 2))),
+    "Expected 'model$domain' to have 1 row, but it had 2",
+    fixed = TRUE)
+  expect_error(
+    mcstate_model(list(density = identity,
+                       parameters = c("a", "b", "c"),
+                       domain = matrix(1:4, 2, 2))),
+    "Expected 'model$domain' to have 3 rows, but it had 2",
+    fixed = TRUE)
+  expect_error(
+    mcstate_model(list(density = identity,
+                       parameters = "a",
+                       domain = matrix(1:4, 1, 4))),
+    "Expected 'model$domain' to have 2 columns, but it had 4",
+    fixed = TRUE)
+})
+
+
+test_that("can use properties to guarantee that a property exists", {
+  m <- list(density = function(x) dnorm(x, log = TRUE), parameters = "a")
+  expect_error(
+    mcstate_model(m, mcstate_model_properties(has_gradient = TRUE)),
+    "Did not find a function 'gradient' within your model")
+  expect_error(
+    mcstate_model(m, mcstate_model_properties(has_direct_sample = TRUE)),
+    "Did not find a function 'direct_sample' within your model")
+  expect_no_error(
+    mcstate_model(m, mcstate_model_properties(has_gradient = FALSE,
+                                              has_direct_sample = FALSE)))
+})
+
+
+test_that("can use properties to ignore properties in a model", {
+  m <- mcstate_model(
+    list(density = identity, parameters = "a", gradient = TRUE),
+    mcstate_model_properties(has_gradient = FALSE))
+  expect_false(m$properties$has_gradient)
+  expect_null(m$gradient)
+
+  m <- mcstate_model(
+    list(density = identity, parameters = "a", direct_sample = TRUE),
+    mcstate_model_properties(has_direct_sample = FALSE))
+  expect_false(m$properties$has_direct_sample)
+  expect_null(m$direct_sample)
+})
+
+
+test_that("require properties are correct type", {
+  expect_error(
+    mcstate_model(
+      list(density = identity, parameters = "a", gradient = TRUE),
+      list(has_gradient = FALSE)),
+    "Expected 'properties' to be a 'mcstate_model_properties' object")
 })
diff --git a/tests/testthat/test-sample.R b/tests/testthat/test-sample.R
index 37d3aaaf..1ad49359 100644
--- a/tests/testthat/test-sample.R
+++ b/tests/testthat/test-sample.R
@@ -132,3 +132,17 @@ test_that("can run more than one chain, in parallel", {
 
   expect_identical(res1, res2)
 })
+
+
+test_that("need a direct sample function in order to start sampling", {
+  model1 <- ex_simple_gamma1()
+  model2 <- mcstate_model(list(density = model1$density,
+                               parameters = model1$parameters,
+                               domain = model1$domain))
+  sampler <- mcstate_sampler_random_walk(vcv = diag(1) * 0.01)
+  expect_no_error(
+    mcstate_sample(model1, sampler, 10))
+  expect_error(
+    mcstate_sample(model2, sampler, 10),
+    "'initial' must be provided with this model")
+})