breedsim-mcp

Breeding-scheme simulation over MCP — returns distributions, never a single stochastic run.

Drives AlphaSimR so an agent can ask what a selection programme would actually gain, with one structural rule: a single simulation run is not a result, and this API will not return one.

Measured on AlphaSimR 2.1.0 — five seeds of an identical three-cycle programme gave mean genetic gain [1.151, 1.841, 1.424, 1.429, 1.473]: sd 0.247 on the very number being reported. Quoting one run to three decimals reports noise with the authority of a measurement. So run_program enforces a replicate floor and returns per-cycle mean, sd and confidence interval. There is no flag that collapses it to a point estimate.

Unofficial. Not affiliated with, endorsed by, or sponsored by the AlphaSimR authors, the University of Edinburgh, or the R Foundation. See NOTICE.

Status

On PyPI — the badge above is the released version, so it cannot go stale the way a number typed here would. Genomic selection included. 5 tools, 57 tests against real AlphaSimR, and 20 mutation checks all confirmed red (docs/MUTATION-CHECKS.md). The simulation is also checked against the breeder's equation R = h²S rather than only against itself (docs/EVAL.md). CI installs R and compiles AlphaSimR, so the suite runs against the real engine on Python 3.11, 3.12 and 3.13 — not against a mock.

Requires mcp 2.x.

Install tax — read this first

Heavier than uv pip install, and the reasons are not negotiable:

• R ≥ 4.3 with a shared library (libR.so)
• AlphaSimR — an Rcpp/RcppArmadillo compile, minutes not seconds
• libtirpc-dev — rpy2 fails to link without it (cannot find -ltirpc)
• rpy2 pinned <3.6 — 3.6 binds R_getVar, which needs R ≥ 4.4

sudo apt-get install -y r-base r-base-dev libtirpc-dev
R -e 'install.packages("AlphaSimR", repos="https://cloud.r-project.org")'
uv add breedsim-mcp

If you build against a conda Python, rpy2 will fail to load libR.so with GLIBCXX_3.4.30 not found — conda ships an older libstdc++ than system libicuuc requires. Use a system or uv-managed interpreter.

Configure your MCP client

The server speaks stdio; the installed console script is breedsim-mcp.

Claude Code

claude mcp add breedsim -- breedsim-mcp

Claude Desktop — add to claude_desktop_config.json:

{
  "mcpServers": {
    "breedsim": {
      "command": "breedsim-mcp"
    }
  }
}

If the executable is not on your PATH, or AlphaSimR lives in a user library, invoke it through uv and pass the library path:

{
  "mcpServers": {
    "breedsim": {
      "command": "uv",
      "args": ["run", "--directory", "/path/to/breedsim-mcp", "breedsim-mcp"],
      "env": { "R_LIBS_USER": "/home/you/R/library" }
    }
  }
}

Verify with list_methods(), which reports the engine versions and whether this process can currently produce reproducible results.

Tools

Typical loop: found_population → run_program → read the CI and the warnings. Comparing two schemes: found_population → compare_programs → read difference. Both run tools take selection_method="phenotypic" or "genomic".

Species

species applies only to runMacs, which carries demographic histories for exactly four: GENERIC, CATTLE, WHEAT, MAIZE (read out of body(runMacs), not the docs). Anything else is refused here rather than failing inside R. Casing does not matter — AlphaSimR upper-cases it, so this does too.

Note the scope that implies: two plants and an animal. Despite the default of MAIZE, this is not a plant-only simulator.

Limits

Every size parameter has a ceiling, reported by list_methods() under limits so a caller can size a request rather than discover the bound by being refused. R runs as a single interpreter here and tool calls are serialised, so one oversized call blocks every other call until it finishes — there is no second worker. The caps are set where a call stops being slow and starts being an outage. For genuinely large jobs, drive AlphaSimR directly rather than through this server.

What run_program returns

Verbatim, for cycles=2, replicates=10, abridged to one cycle:

{
  "session_id": "bs-...",
  "replicates": 10,
  "cycles": [
    {
      "cycle": 2,
      "genetic_gain": {
        "mean": 1.5941703785773211,
        "sd": 0.1414620049281876,
        "ci_low": 1.4929815869737013,
        "ci_high": 1.695359170180941,
        "n": 10
      },
      "genetic_variance": {
        "mean": 0.5149269761002959,
        "sd": 0.15228915685789685,
        "ci_low": 0.4059934446929936,
        "ci_high": 0.6238605075075982,
        "n": 10
      }
    }
  ],
  "reproducible": true,
  "recipe": {
    "generator": "quickHaplo",
    "seed": 1,
    "n_select": 10,
    "n_cross": 60,
    "base_seed": 1000
  },
  "warnings": []
}

There is no value field anywhere. Intervals use t critical values rather than a normal 1.96, because at n = 5–10 the normal understates the interval — the wrong direction to be wrong in when the interval exists to be honest.

Comparing two programmes

Do not call run_program twice and compare the means. Use compare_programs, which pairs the two arms on the same seeds — replicate i of A and replicate i of B start from identical founders under an identical seed — and differences them within each pair, so the shared luck of that seed cancels instead of being counted twice.

Read difference and favours. favours is null when the interval contains zero, which means the two programmes are not distinguishable at that replicate count; the larger mean is then not the better programme.

Here is why the pairing earns its keep. Verbatim, selecting 12 of 100 against 18 of 100, final cycle of two, ten replicates:

{
  "programs": {
    "a": { "label": "A", "n_select": 12, "n_cross": 100 },
    "b": { "label": "B", "n_select": 18, "n_cross": 100 }
  },
  "cycles": [
    {
      "cycle": 2,
      "a_genetic_gain": {
        "mean": 2.046227841067686,
        "ci_low": 1.9001469542399823,
        "ci_high": 2.1923087278953903,
        "n": 10
      },
      "b_genetic_gain": {
        "mean": 1.730473406465538,
        "ci_low": 1.5520780571733739,
        "ci_high": 1.9088687557577022,
        "n": 10
      },
      "difference": {
        "mean": 0.31575443460214814,
        "ci_low": 0.10028439648886733,
        "ci_high": 0.5312244727154289,
        "n": 10
      }
    }
  ],
  "favours": "a",
  "intervals_overlap": true,
  "warnings": [{ "code": "overlap_but_different", "message": "..." }]
}

The two per-programme intervals overlap — A spans 1.900–2.192, B spans 1.552–1.909 — so reading them side by side says "no difference". The paired difference says otherwise: [+0.100, +0.531], entirely above zero. Pairing cancels the seed-to-seed variation that made both individual intervals wide, so it resolves a contrast that eyeballing the overlap cannot. That is what overlap_but_different is for.

Two overlapping confidence intervals do not imply no difference. This is the single easiest way to get a breeding comparison wrong, and it is why the tool reports a difference rather than two numbers.

Genomic selection — and the trap under it

selection_method="genomic" fits RRBLUP to the marker genotypes each cycle and selects on the estimated breeding value instead of the phenotype. It needs a SNP chip, which is a founding decision:

found_population(generator="runMacs", n_snp_per_chr=50)  # note the generator
run_program(session_id, selection_method="genomic")

Note the generator, because this is where genomic selection goes quietly wrong. Markers predict a trait only through linkage disequilibrium with the causal loci — that is the whole mechanism. And quickHaplo, the default generator and the only reproducible one, has none:

Adjacent markers in quickHaplo are no more correlated than randomly chosen distant ones — it samples haplotypes with no coalescent history, so there is no linkage to learn from. Every found_population call with a chip therefore returns a measured linkage_disequilibrium block, and a no_linkage_disequilibrium warning when the ratio says the markers are uninformative.

So on this engine, reproducibility and genomic realism cannot be had at the same time. quickHaplo reproduces and cannot support genomic selection; runMacs supports it and does not reproduce in a long-lived process. That is a real constraint of the underlying simulator, and the server states it rather than letting you find it as a wrong answer.

Why the guard measures LD instead of accuracy

The obvious alternative — fit the model, warn if accuracy is poor — cannot do the job. Measured at 20 replicates, 200 individuals, three cycles:

A population with no linkage disequilibrium at all reaches 0.208, and at 10% selection it beats the population that has real LD. Accuracy also climbs every cycle in both. Neither is a paradox: out-of-sample accuracy in a closed population conflates LD with the causal loci and plain relatedness between training and target individuals. As descendants of a few selected parents fill the population, markers predict by tracking pedigree — and quickHaplo's mutually uncorrelated markers tag pedigree more efficiently than runMacs' markers, which are partly redundant with each other precisely because they are in LD.

An accuracy threshold would therefore wave through the exact population it claimed to catch. Only the LD measurement discriminates, so that is what gates the warning. Accuracy is still reported on every genomic cycle, measured out-of-sample on progeny the model never saw — reporting its in-sample fit instead would have read 0.448 where the truth was 0.097. Read it as a property of the model in front of you, not as proof that genomic selection is working for the reason you assume.

Warnings

None of these withhold results. Outputs are always distributions, so you can already see when an answer is too noisy to use — they explain rather than refuse.

Reproducibility

Two independent things break it, and both were measured:

The runMacs case is subtler than "it ignores the seed". Its MaCS RNG is seeded once per R session and advances across calls, so the same call sequence reproduces in a fresh process while a repeat call inside one session does not. Since this server is a long-lived process, the repeat case is the one you hit — hence reproducible: false.

The server pins threads at import, before rpy2 loads, and reports whether it succeeded. With founders fixed and threads pinned, seed 7 reproduces exactly: meanG=2.01451853, twice.

Limitations

Multi-trait selection is supported: pass h2 as a list to build several traits and index_weights to select on a weighted index. compare_programs remains single-trait — one paired verdict needs one criterion, and with several traits that criterion is the index, which AlphaSimR does not report a gain for.

No G×E, no optimal contribution selection, no crossing-block optimisation, no genotype-matrix export. Genomic selection is RRBLUP only — the RRBLUP_D, _GCA and _SCA variants are not exposed, because they estimate dominance and combining-ability effects that an additive-only addTraitA architecture does not contain.

The prediction_accuracy_low and no_linkage_disequilibrium warnings are advisory: like every other warning here, they explain rather than refuse.

Sessions are in-memory and capped (8, LRU); they do not survive a restart.

Licence

GPL-3.0-or-later — because this imports rpy2, which is GPLv2+. AlphaSimR itself is MIT; the licence here is about rpy2, not about AlphaSimR. See NOTICE.

More

• CHANGELOG.md — what changed, and why
• docs/MUTATION-CHECKS.md — every guard disabled on purpose, and the test that went red for it