DAFJS09 |  | [28,845, 42,537†] 32.19% | 130,367 | 100 | 3,600.0 |  | [16,016, 36,539] 56.17% | 1,722,093 | 3,600.0
DAFJS10 |  | [30,021, 46,257†] 35.1% | 63,949 | 10 | 3,600.0 |  | [16,354, 40,354] 59.47% | 1,779,988 | 3,600.0
DAFJS11 |  | [29,853, 56,847†] 47.49% | 55,377 | 0 | 3,600.0 |  | [29,499, 55,646] 46.99% | 205,807 | 3,600.0
DAFJS12 |  | [26,886, 55,575†] 51.62% | 99,798 | 0 | 3,600.0 |  | [23,461, 55,575†] 57.78% | 104,236 | 3,600.0
DAFJS13 |  | [39,200, 51,317†] 23.61% | 50,909 | 66 | 3,600.0 |  | [14,433, 48,176] 70.04% | 880,264 | 3,600.0
DAFJS14 |  | [43,922, 62,724†] 29.98% | 36,836 | 13 | 3,600.0 |  | [16,518, 54,861] 69.89% | 769,749 | 3,600.0
DAFJS15 |  | [0, 54,353†] 100% | 55,931 | 0 | 3,600.0 |  | [22,324, 53,549] 58.31% | 97,946 | 3,600.0
DAFJS16 |  | [28,823, 65,550†] 56.03% | 92,513 | 0 | 3,600.0 |  | [28,167, 65,257] 56.84% | 135,220 | 3,600.0
DAFJS17 |  | [42,492, 63,177†] 32.74% | 33,771 | 0 | 3,600.0 |  | [13,683, 56,330] 75.71% | 447,082 | 3,600.0
DAFJS18 |  | [45,630, 61,831†] 26.2% | 29,131 | 23 | 3,600.0 |  | [14,813, 59,787] 75.22% | 1,007,448 | 3,600.0
DAFJS19 |  | [27,728, 58,001†] 52.19% | 47,626 | 17 | 3,600.0 |  | [24,250, 44,486] 45.49% | 881,311 | 3,600.0
DAFJS20 |  | [37,194, 63,588†] 41.51% | 42,948 | 0 | 3,600.0 |  | [19,568, 63,065] 68.97% | 332,996 | 3,600.0
DAFJS21 |  | [36,582, 63,202†] 42.12% | 54,908 | 0 | 3,600.0 |  | [22,142, 62,438] 64.54% | 179,594 | 3,600.0
DAFJS22 |  | [31,266, 56,115†] 44.28% | 72,153 | 0 | 3,600.0 |  | [19,800, 56,115†] 64.72% | 146,289 | 3,600.0
DAFJS23 |  | [23,845, 47,616†] 49.92% | 40,481 | 0 | 3,600.0 |  | [21,422, 42,266] 49.32% | 488,429 | 3,600.0
DAFJS24 |  | [26,527, 49,019†] 45.88% | 56,184 | 0 | 3,600.0 |  | [21,763, 48,525] 55.15% | 346,533 | 3,600.0
DAFJS25 |  | [0, 66,262†] 100% | 111,844 | 0 | 3,600.0 |  | [24,843, 66,200] 62.47% | 88,055 | 3,600.0
DAFJS26 |  | [0, 72,226†] 100% | 77,334 | 0 | 3,600.0 |  | [24,182, 71,874] 66.36% | 101,029 | 3,600.0
DAFJS27 |  | [0, 64,797†] 100% | 64,933 | 0 | 3,600.0 |  | [21,584, 64,781] 66.68% | 145,287 | 3,600.0
DAFJS28 |  | [25,167, 52,639†] 52.19% | 58,308 | 0 | 3,600.0 |  | [24,712, 50,995] 51.54% | 193,331 | 3,600.0
DAFJS29 |  | [29,463, 59,792†] 50.72% | 66,255 | 0 | 3,600.0 |  | [27,497, 57,222] 51.95% | 246,809 | 3,600.0
DAFJS30 |  | [27,084, 55,015†] 50.77% | 91,519 | 0 | 3,600.0 |  | [21,222, 54,467] 61.04% | 306,321 | 3,600.0
YFJS01 |  | [47,837, 84,152†] 43.15% | 768,698 | 1,077 | 3,600.0 |  | [39,902, 63,486] 37.15% | 4,653,807 | 3,600.0
YFJS02 |  | [59,193, 66,853†] 11.46% | 1,587,076 | 2,455 | 3,600.0 |  | [46,660, 66,853†] 30.21% | 16,559,119 | 3,600.0
YFJS03 |  | 30,073 | 2,666,461 | 8,869 | 651.0 |  | 30,073 | 1,859,541 | 89.8
YFJS04 |  | [29,362, 35,810] 18.01% | 8,026,009 | 27,106 | 3,600.0 |  | [21,857, 32,670] 33.1% | 44,018,889 | 3,600.0
YFJS05 |  | [34,400, 40,221] 14.47% | 6,938,357 | 20,891 | 3,600.0 |  | [24,921, 37,682] 33.86% | 37,757,176 | 3,600.0
YFJS06 |  | [31,101, 38,684] 19.6% | 3,339,460 | 7,174 | 3,600.0 |  | [19,541, 36,812] 46.92% | 17,596,410 | 3,600.0
YFJS07 |  | [31,513, 42,140] 25.22% | 2,383,275 | 7,262 | 3,600.0 |  | [26,063, 37,664] 30.8% | 11,920,798 | 3,600.0
YFJS08 |  | 30,184 | 8,684,046 | 11,067 | 2,401.4 |  | 30,184 | 15,779,708 | 1,242.7
YFJS09 |  | [14,785, 23,430] 36.9% | 330,752 | 106 | 3,600.0 |  | [11,324, 21,796] 48.05% | 5,536,683 | 3,600.0
YFJS10 |  | [30,116, 37,720] 20.16% | 8,688,281 | 15,350 | 3,600.0 |  | [22,264, 34,624] 35.7% | 38,762,863 | 3,600.0
YFJS11 |  | [35,760, 51,469†] 30.52% | 653,855 | 652 | 3,600.0 |  | [30,073, 44,294] 32.11% | 16,060,438 | 3,600.0
YFJS12 |  | [32,446, 59,898†] 45.83% | 795,705 | 1,181 | 3,600.0 |  | [26,939, 43,790] 38.48% | 13,545,313 | 3,600.0
YFJS13 |  | [25,469, 45,084†] 43.51% | 466,368 | 549 | 3,600.0 |  | [21,361, 35,266] 39.43% | 8,512,041 | 3,600.0
YFJS14 |  | [68,791, 109,469†] 37.16% | 74,831 | 0 | 3,600.0 |  | [68,423, 98,323] 30.41% | 1,965,536 | 3,600.0
YFJS15 |  | [65,453, 112,424†] 41.78% | 78,095 | 0 | 3,600.0 |  | [64,923, 105,088] 38.22% | 5,880,443 | 3,600.0
YFJS16 |  | [63,661, 106,855†] 40.42% | 56,452 | 0 | 3,600.0 |  | [62,865, 103,358] 39.18% | 1,474,916 | 3,600.0
YFJS17 |  | [0, 85,736†] 100% | 84,708 | 0 | 3,600.0 |  | [51,540, 85,611] 39.8% | 277,616 | 3,600.0
YFJS18 |  | [0, 99,429†] 100% | 85,143 | 0 | 3,600.0 |  | [55,101, 97,691] 43.6% | 256,080 | 3,600.0
YFJS19 |  | [0, 89,561†] 100% | 109,068 | 0 | 3,600.0 |  | [41,690, 89,561†] 53.45% | 205,157 | 3,600.0
YFJS20 |  | [0, 91,958†] 100% | 111,180 | 0 | 3,600.0 |  | [43,711, 91,958†] 52.47% | 127,010 | 3,600.0
Mean |  | [28348.02, 57881.72] 51.02% | 1,013,398.90 | 2,224.32 | 3,517.05 |  | [27680.84, 53967.68] 48.71% | 5,897,637.44 | 3,482.65

Table 15: Solutions found and computational cost of solving the large-sized
instances with learning rate $\alpha=0.2$ using CPLEX and CP Optimizer with
warm start.

instance |  | IBM ILOG CPLEX |  | IBM ILOG CP Optimizer
---|---|---|---|---
| makespan | Effort measurement |  | makespan | Effort measurement
| #iterations | #B&B Nodes | CPU |  | #branches | CPU
DAFJS01 |  | [13,085, 22,616†] 42.14% | 1,117,444 | 1,730 | 3,600.0 |  | [10,730, 20,247] 47% | 7,775,777 | 3,600.0
DAFJS02 |  | [16,801, 26,183†] 35.83% | 1,515,028 | 3,091 | 3,600.0 |  | [10,500, 22,147] 52.59% | 18,313,536 | 3,600.0
DAFJS03 |  | [22,719, 43,964†] 48.32% | 87,465 | 6 | 3,600.0 |  | [21,152, 42,677] 50.44% | 653,668 | 3,600.0
DAFJS04 |  | [24,927, 43,281†] 42.41% | 150,187 | 31 | 3,600.0 |  | [23,278, 40,023] 41.84% | 2,127,914 | 3,600.0
DAFJS05 |  | [20,474, 38,237†] 46.46% | 298,694 | 514 | 3,600.0 |  | [12,401, 29,897] 58.52% | 3,366,879 | 3,600.0
DAFJS06 |  | [20,171, 38,007†] 46.93% | 103,184 | 100 | 3,600.0 |  | [12,004, 28,181] 57.4% | 1,790,536 | 3,600.0
DAFJS07 |  | [16,887, 41,285†] 59.1% | 57,392 | 0 | 3,600.0 |  | [15,801, 40,857] 61.33% | 311,704 | 3,600.0
DAFJS08 |  | [21,334, 46,323†] 53.95% | 84,392 | 0 | 3,600.0 |  | [20,629, 42,452] 51.41% | 289,066 | 3,600.0
DAFJS09 |  | [22,185, 41,418†] 46.44% | 176,125 | 14 | 3,600.0 |  | [11,423, 33,332] 65.73% | 1,702,714 | 3,600.0
DAFJS10 |  | [22,463, 38,771†] 42.06% | 71,903 | 14 | 3,600.0 |  | [11,408, 32,284] 64.66% | 1,290,793 | 3,600.0
DAFJS11 |  | [20,011, 47,805†] 58.14% | 81,360 | 0 | 3,600.0 |  | [19,709, 47,731] 58.71% | 193,904 | 3,600.0
DAFJS12 |  | [18,895, 45,735†] 58.69% | 76,070 | 0 | 3,600.0 |  | [15,608, 44,045] 64.56% | 196,960 | 3,600.0
DAFJS13 |  | [29,210, 42,319†] 30.98% | 65,284 | 26 | 3,600.0 |  | [9,929, 40,101] 75.24% | 1,351,794 | 3,600.0
DAFJS14 |  | [32,979, 49,459†] 33.32% | 47,535 | 37 | 3,600.0 |  | [11,318, 45,078] 74.89% | 763,657 | 3,600.0
DAFJS15 |  | [0, 49,206†] 100% | 71,721 | 0 | 3,600.0 |  | [14,671, 48,366] 69.67% | 131,171 | 3,600.0
DAFJS16 |  | [20,218, 56,727†] 64.36% | 97,311 | 0 | 3,600.0 |  | [18,661, 56,727†] 67.1% | 171,678 | 3,600.0
DAFJS17 |  | [29,862, 54,527†] 45.23% | 29,116 | 0 | 3,600.0 |  | [9,240, 48,773] 81.06% | 499,025 | 3,600.0
DAFJS18 |  | [32,271, 51,823†] 37.73% | 29,115 | 12 | 3,600.0 |  | [10,047, 47,845] 79% | 670,899 | 3,600.0
DAFJS19 |  | [21,074, 46,880†] 55.05% | 44,072 | 9 | 3,600.0 |  | [16,675, 43,407] 61.58% | 654,923 | 3,600.0
DAFJS20 |  | [25,682, 51,686†] 50.31% | 58,048 | 0 | 3,600.0 |  | [13,124, 51,061] 74.3% | 333,744 | 3,600.0
DAFJS21 |  | [0, 53,010†] 100% | 66,397 | 0 | 3,600.0 |  | [14,664, 50,104] 70.73% | 234,319 | 3,600.0
DAFJS22 |  | [21,312, 45,005†] 52.65% | 85,113 | 0 | 3,600.0 |  | [12,934, 45,005†] 71.26% | 146,342 | 3,600.0
DAFJS23 |  | [18,008, 39,249†] 54.12% | 41,167 | 0 | 3,600.0 |  | [14,782, 39,145] 62.24% | 652,256 | 3,600.0
DAFJS24 |  | [18,974, 44,851†] 57.7% | 61,919 | 0 | 3,600.0 |  | [14,688, 44,462] 66.97% | 255,869 | 3,600.0
DAFJS25 |  | [0, 51,964†] 100% | 71,177 | 0 | 3,600.0 |  | [16,174, 51,964†] 68.87% | 83,334 | 3,600.0
DAFJS26 |  | [0, 57,182†] 100% | 103,993 | 0 | 3,600.0 |  | [15,820, 57,182†] 72.33% | 123,884 | 3,600.0
DAFJS27 |  | [25,010, 58,138†] 56.98% | 107,828 | 0 | 3,600.0 |  | [14,138, 57,409] 75.37% | 139,297 | 3,600.0
DAFJS28 |  | [18,433, 42,898†] 57.03% | 63,459 | 0 | 3,600.0 |  | [16,781, 40,331] 58.39% | 611,659 | 3,600.0
DAFJS29 |  | [22,832, 51,606†] 55.76% | 58,895 | 0 | 3,600.0 |  | [18,452, 51,156] 63.93% | 194,564 | 3,600.0
DAFJS30 |  | [17,829, 43,859†] 59.35% | 82,952 | 0 | 3,600.0 |  | [14,275, 42,315] 66.26% | 281,586 | 3,600.0
YFJS01 |  | [39,168, 66,116] 40.76% | 560,775 | 1,035 | 3,600.0 |  | [29,808, 58,170] 48.76% | 6,636,736 | 3,600.0
YFJS02 |  | [48,544, 59,211] 18.02% | 1,188,834 | 2,163 | 3,600.0 |  | [35,526, 59,001] 39.79% | 12,010,945 | 3,600.0
YFJS03 |  | 27,686 | 1,406,947 | 4,582 | 337.1 |  | 27,686 | 1,278,286 | 66.9
YFJS04 |  | [25,715, 32,191] 20.12% | 8,529,348 | 26,006 | 3,600.0 |  | [18,161, 29,692] 38.84% | 42,209,341 | 3,600.0
YFJS05 |  | [31,090, 34,905] 10.93% | 8,010,184 | 28,009 | 3,600.0 |  | [19,314, 35,629] 45.79% | 40,473,041 | 3,600.0
YFJS06 |  | [26,650, 39,813] 33.06% | 2,057,905 | 4,922 | 3,600.0 |  | [14,714, 33,822] 56.5% | 14,236,397 | 3,600.0
YFJS07 |  | [27,041, 38,134] 29.09% | 3,907,838 | 10,102 | 3,600.0 |  | [20,397, 33,059] 38.3% | 12,005,380 | 3,600.0
YFJS08 |  | [27,020, 28,638] 5.65% | 11,006,083 | 13,663 | 3,600.0 |  | [17,787, 28,192] 36.91% | 41,659,118 | 3,600.0
YFJS09 |  | [12,166, 21,816] 44.23% | 439,398 | 350 | 3,600.0 |  | [9,525, 20,457] 53.44% | 4,119,430 | 3,600.0
YFJS10 |  | [28,406, 33,622] 15.51% | 8,543,764 | 16,090 | 3,600.0 |  | [17,808, 32,159] 44.63% | 30,351,408 | 3,600.0
YFJS11 |  | [30,906, 45,051†] 31.4% | 770,066 | 1,119 | 3,600.0 |  | [23,170, 41,095] 43.62% | 7,467,430 | 3,600.0
YFJS12 |  | [26,035, 48,493] 46.31% | 381,810 | 116 | 3,600.0 |  | [21,335, 39,039] 45.35% | 7,124,912 | 3,600.0
YFJS13 |  | [23,227, 40,119†] 42.1% | 812,817 | 1,074 | 3,600.0 |  | [16,668, 30,711] 45.73% | 4,335,509 | 3,600.0
YFJS14 |  | [49,960, 92,457†] 45.96% | 83,456 | 0 | 3,600.0 |  | [49,252, 91,042] 45.9% | 1,297,229 | 3,600.0
YFJS15 |  | [47,574, 96,547†] 50.72% | 77,782 | 0 | 3,600.0 |  | [46,947, 89,100] 47.31% | 824,181 | 3,600.0
YFJS16 |  | [47,654, 92,811†] 48.65% | 55,873 | 0 | 3,600.0 |  | [45,648, 87,080] 47.58% | 1,706,868 | 3,600.0
YFJS17 |  | [0, 73,682†] 100% | 65,088 | 0 | 3,600.0 |  | [34,763, 73,682†] 52.82% | 110,226 | 3,600.0
YFJS18 |  | [0, 87,059†] 100% | 72,375 | 0 | 3,600.0 |  | [37,042, 87,059†] 57.45% | 111,311 | 3,600.0
YFJS19 |  | [0, 74,431†] 100% | 62,604 | 0 | 3,600.0 |  | [27,976, 74,431†] 62.41% | 138,828 | 3,600.0
YFJS20 |  | [0, 72,481†] 100% | 133,020 | 0 | 3,600.0 |  | [29,376, 72,250] 59.34% | 138,581 | 3,600.0
Mean |  | [21849.76, 49385.54] 51.47% | 1,063,406.26 | 2,296.30 | 3,534.74 |  | [19678.38, 46553.2] 57.73% | 5,470,972.18 | 3,529.34

Table 16: Solutions found and computational cost of solving the large-sized
instances with learning rate $\alpha=0.3$ using CPLEX and CP Optimizer with
warm start.
---
Figure 6: Graphical representation of the information contained in Tables
6–16. The figures at the top, middle and bottom correspond to $\alpha$ equal
to $0.1$, $0.2$, and $0.3$, respectively. Each figure shows, for each large-
sized size instance, the solution found by solving the MILP and CP models with
and without warm start, as well as the lower bounds found in these 4 cases.
The numbers from 1 to 50 on the abscissa axis correspond to the 50 large-sized
instances in the order presented in all the tables. That is, the instances
from 1 to 30 correspond to the instances DAFJS01, $\dots$, DAFJS30 and the
instances from 31 to 50 correspond to the instances YFJS01, $\dots$, YFJS20.

##  6 Conclusions

In this work, we addressed the FJS scheduling problem with sequencing
flexibility and position-based learning effect. To the authors’ knowledge,
this combination, with potential application in a wide range of real-world
industrial environments, has never been addressed in the literature. As a
first step towards its efficient and effective solution, we introduced a set
of 110 instances that transform into 330 instances by varying the learning
rate $\alpha\in\\{0.1,0.2,0.3\\}$. By introducing MILP and CP models, an exact
solver aided by constructive heuristics was able to provide 183 proven optimal
solutions. Instances, their solutions, models and constructive heuristics are
all available for download at https://github.com/kennedy94/FJS. We expect this
benchmark test-set to guide the introduction and evaluation of new effective
heuristic and metaheuristic approaches for the considered problem. In fact,
this is the current line of research of the authors.

Conflict of interest statement: On behalf of all authors, the corresponding
author states that there is no conflict of interest.

Data availability: The datasets generated during and/or analysed during the
current study are available in the GitHub repository,
https://github.com/kennedy94/FJS.

## References

  * [1] R. Alvarez-Valdes, A. Fuertes, J. M. Tamarit, G. Giménez, and R. Ramos.  A heuristic to schedule flexible job-shop in a glass factory.  European Journal of Operational Research, 165(2):525–534, 2005.
  * [2] J. L. Andrade-Pineda, D. Canca, P. L. Gonzalez-R, and M. Calle.  Scheduling a dual-resource flexible job shop with makespan and due date-related criteria.  Annals of Operations Research, 291(1):5–35, 2020.
  * [3] A. Azzouz, M. Ennigrou, and L. Ben Said.  Scheduling problems under learning effects: classification and cartography.  International Journal of Production Research, 56(4):1642–1661, 2017.
  * [4] E. G. Birgin, P. Feofiloff, C. G. Fernandes, E. L. De Melo, M. T. I. Oshiro, and D. P. Ronconi.  A MILP model for an extended version of the flexible job shop problem.  Optimization Letters, 8(4):1417–1431, 2014.
  * [5] E. G. Birgin, J. E. Ferreira, and D. P. Ronconi.  List scheduling and beam search methods for the flexible job shop scheduling problem with sequencing flexibility.  European Journal of Operational Research, 247(2):421–440, 2015.
  * [6] D. Biskup.  Single-machine scheduling with learning considerations.  European Journal of Operational Research, 115(1):173–178, 1999.
  * [7] D. Biskup.  A state-of-the-art review on scheduling with learning effects.  European Journal of Operational Research, 188(2):315–329, 2008.
  * [8] Z.-C. Cao, C.-R. Lin, and M.-C. Zhou.  A knowledge-based cuckoo search algorithm to schedule a flexible job shop with sequencing flexibility.  IEEE Transactions on Automation Science and Engineering, 18(1):56–69, 2021.
  * [9] T. C. E. Cheng and G. Wang.  Single machine scheduling with learning effect considerations.  Annals of Operations Research, 98(1/4):273–290, 2000.
  * [10] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein.  Introduction to Algorithms.  The MIT Press, Cambridge, MA, USA, 4th edition, 2022.
  * [11] P. Y. Gan and K. S. Lee.  Scheduling of flexible-sequenced process plans in a mould manufacturing shop.  The International Journal of Advanced Manufacturing Technology, 20(3):214–222, 2002.
  * [12] J. Gao, M. Gen, and L. Sun.  Scheduling jobs and maintenances in flexible job shop with a hybrid genetic algorithm.  Journal of Intelligent Manufacturing, 17(4):493–507, 2006.
  * [13] M. R. Garey, D. S. Johnson, and R. Sethi.  The complexity of flowshop and jobshop scheduling.  Mathematics of Operations Research, 1(2):117–129, 1976.
  * [14] J. N. D. Gupta and S. K. Gupta.  Single facility scheduling with nonlinear processing times.  Computers and Industrial Engineering, 14(4):387–393, 1988.
  * [15] A. Janiak, T. Krysiak, and R. Trela.  Scheduling problems with learning and ageing effects: A survey.  Decision Making in Manufacturing and Services, 5(1):19–36, 2011.
  * [16] G. A. Kasapidis, S. Dauzère‐Pérès, D. C. Paraskevopoulos, P. P. Repoussis, and C. D. Tarantilis.  On the multiresource flexible job‐shop scheduling problem with arbitrary precedence graphs.  Production and Operations Management, 32(7):2322–2330, 2023.
  * [17] G. A. Kasapidis, D. C. Paraskevopoulos, P. P. Repoussis, and C. D. Tarantilis.  Flexible job shop scheduling problems with arbitrary precedence graphs.  Production and Operations Management, 30(11):4044–4068, 2021.
  * [18] Y. K. Kim, K. Park, and J. Ko.  A symbiotic evolutionary algorithm for the integration of process planning and job shop scheduling.  Computers and Operations Research, 30(8):1151–1171, 2003.
  * [19] P. Laborie, J. Rogerie, P. Shaw, and P. Vilím.  IBM ILOG CP optimizer for scheduling: 20+ years of scheduling with constraints at ibm/ilog.  Constraints, 23(2):210–250, 2018.
  * [20] J. Y.-T. Leung, H. Li, and M. Pinedo.  Order scheduling in an environment with dedicated resources in parallel.  Journal of Scheduling, 8(5):355–386, 2005.
  * [21] W. T. Lunardi, E. G. Birgin, P. Laborie, D. P. Ronconi, and H. Voos.  Mixed integer linear programming and constraint programming models for the online printing shop scheduling problem.  Computers and Operations Research, 123:105020, 2020.
  * [22] W. T. Lunardi, E. G. Birgin, D. P. Ronconi, and H. Voos.  Metaheuristics for the online printing shop scheduling problem.  European Journal of Operational Research, 293(2):419–441, 2021.
  * [23] G. Vilcot and J.-C. Billaut.  A tabu search and a genetic algorithm for solving a bicriteria general job shop scheduling problem.  European Journal of Operational Research, 190(2):398–411, 2008.
  * [24] A. Vital-Soto, A. Azab, and M. F. Baki.  Mathematical modeling and a hybridized bacterial foraging optimization algorithm for the flexible job-shop scheduling problem with sequencing flexibility.  Journal of Manufacturing Systems, 54:74–93, 2020.
  * [25] H. M. Wagner.  An integer linear-programming model for machine scheduling.  Naval Research Logistics Quarterly, 6(2):131–140, 1959.
  * [26] J. M. Wilson.  Alternative formulations of a flow-shop scheduling problem.  Journal of the Operational Research Society, 40(4):395–399, 1989.
  * [27] L. Yu, C. Zhu, J. Shi, and W. Zhang.  An extended flexible job shop scheduling model for flight deck scheduling with priority, parallel operations, and sequence flexibility.  Scientific Programming, 2017:1–15, 2017.


# Multi-dimensional Energy Limitation in Sphere Shaping for Nonlinear
Interference Noise Mitigation

Jingtian Liu  Télécom Paris, Institut Polytechnique de Paris
19 place Marguerite Perey, 91120 Palaiseau, FRANCE
<EMAIL_ADDRESS>Élie Awwad  Yves Jaouën  Télécom Paris,
Institut Polytechnique de Paris
19 place Marguerite Perey, 91120 Palaiseau, FRANCE
<EMAIL_ADDRESS>Télécom Paris, Institut Polytechnique de Paris
19 place Marguerite Perey, 91120 Palaiseau, FRANCE
<EMAIL_ADDRESS>
###### Abstract

We propose Four-Dimensional (4D) energy limit enumerative sphere shaping (ESS)
of $M$-QAM signaling to minimize rate loss and improve the transmission
performance over non-linear WDM optical-fiber systems. Simulation results show
that the proposed scheme outperforms the conventional ESS by $0.19$
bit/4D-symbol in achievable information rate over a $205$-km single-span link
and a WDM transmission of five polarization-division-multiplexed channels with
$400$-Gbit/s net rate per channel. We also study the achieved performance over
several shaping block lengths and show that the achieved gains do not scale
well over multi-span systems.

###### Index Terms:

Probabilistic constellation shaping, enumerative sphere shaping, coherent
transmission systems, modulation and coding schemes, non-linear Kerr effect.

##  I Introduction

Probabilistic amplitude shaping (PAS) has demonstrated excellent performance
over linear AWGN channels. PAS can be implemented through several approaches,
among which we cite enumerative sphere shaping (ESS) [1] and constant
composition [2] distribution matching (CCDM). When it comes to dealing with
non-linear interference (NLI) in optical fiber transmission, PAS may be more
susceptible to performance degradation [3]. Recently, several works aimed at
improving the non-linear performance of PAS schemes by limiting the
1-dimensional (1D) amplitude variance or kurtosis of ESS schemes [4, 5] or by
limiting the 2D energy dispersion index (EDI) of CCDM schemes [6]. Inspired by
these advancements, we propose a new distribution matching approach to
alleviate the distortions introduced by Kerr effects in optical fiber systems.
Our design targets lowering 4D energy variations where the four dimensions
refer to the in-phase and quadrature parts of two orthogonal polarization
tributaries of the transmitted signal. The proposed scheme is based on
polarization-division-multiplexed (PDM) quadrature-and-amplituded (QAM)
signaling along with ESS for amplitude shaping. ESS is chosen for its reduced
rate loss when it is applied over short shaping block lengths compared to CCDM
[1].

In this paper, we propose a 4D energy limit over ESS shaping to achieve either
an improved non-linear tolerance or a lower rate loss than state-of-the-art
kurtosis-ESS (K-ESS) [4] and band-limited ESS (B-ESS) [5], and this is done by
tuning the design parameters of the new scheme. Simulation results show that
the proposed shaping outperforms conventional ESS in two different scenarios
when applied over $64$-QAM. We measure a gain of $0.19$ bit per 4D-symbol in
the achievable information rate (AIR) over a $205$-km single-span standard
single-mode-fiber (SSMF) link and a WDM transmission of five channels with a
net rate of $400$ Gbit/s per wavelegnth at $50$ GBaud and a block length of
$N=108$. Smaller gains are observed for larger block lengths over a $5$-span
$400$-km link and a WDM transmission of five $110$ GBaud channels.

##  II Principle

ESS can be generated and applied over $1$, $2$, $4$ or even more dimensions.
In Fig. 3 from [7], the authors show an illustration of 1D, 2D, and 4D symbol
mapping strategies of ESS schemes. For 1D mapping, 4D symbols are shaped using
$4$ independent amplitude sequences (one ESS per dimension). For 2D mapping,
each polarization, namely $X$ and $Y$, uses one ESS-shaped amplitude sequence.
For 4D mapping, a single-shaped sequence is used for a wavelength channel. The
latter mapping offers the advantage of reducing the probability of having
equal amplitudes at the same time slot; hence, high peak power values appear
less frequently and lead to lower non-linear distortions. This advantage was
demonstrated in several works [7, 4, 3].

Figure 1: 4D-Band-limited-trellis using four amplitudes $\left\\{1,3,5,7\right\\}$, $E_{\mathrm{max}}=48$, and block length $N=8$. TABLE I: Metrics for five different ESS schemes with $N=108$, rate = $1.5$ bits/1D amplitude | ESS | 1D-BL-ESS | 4D-BL-ESS (Linear) | 4D-BL-ESS (Nonlinear) | K-ESS
---|---|---|---|---|---
$E_{\mathrm{max}}$ | 860 | 996 | 948 | 996 | 1156
$\overline{E}_{1D}$ | 7.85 | 9.14 | 8.63 | 9.09 | 8.27
$\mathrm{\sigma^{2}_{1D}}$ | 1.75 | 1.49 | 1.52 | 1.45 | 1.16
$\mathrm{\mu_{4,2D}}$ | 1.87 | 1.71 | 1.69 | 1.64 | 1.57
$\mathrm{EDI}~{}(W_{1}=42)$ | 3.13 | 0.65 | 0.72 | 0.48 | 3.34
$\overline{\mu}_{4,E_{4D}}(W_{2}=54)$ | 3.63 | 2.8 | 2.84 | 2.69 | 3.17
$\mathrm{Rate~{}loss~{}[bit/1D]}$ | 0.03 | 0.13 | 0.097 | 0.13 | 0.07

While both K-ESS [4] and B-ESS [5] (that we call 1D-BL-ESS in the following)
use 4D mapping, they introduce energy constraints over a single dimension to
limit the variations of each 1D amplitude. However, the strength of non-linear
effects depends on the variations of energy levels of the 4D symbols in
consecutive time slots. Hence, by constraining energy variations of 4D
symbols, we can improve the performance and achieve a better trade-off between
linear and non-linear gains. The following toy example illustrates this
advantage. Consider an ESS with a block length $N=8$ and two amplitude
sequences $(3,3,3,3,1,1,1,1)$, highlighted in yellow in Fig. 1, and
$(1,3,1,3,1,3,1,3)$, highlighted in green. These sequences have the same 1D
average energy, variance, and kurtosis, and generate the two consecutive
4D-symbols: $[\left(\pm 3\pm 3i,\pm 3\pm 3i\right)^{T},\left(\pm 1\pm 1i,\pm
1\pm 1i\right)^{T}]$ and $[\left(\pm 1\pm 3i,\pm 1\pm 3i\right)^{T},\left(\pm
1\pm 3i,\pm 1\pm 3i\right)^{T}]$ respectively where the first symbol is
attributed to the polarization $X$ and the second symbol to the polarization
$Y$ in each pair, and $(\cdot)^{T}$ is the transpose of $(\cdot)$. Symbols
from the first sequence show a large temporal 4D energy fluctuation from $36$
to $4$ while symbols from the second one have a constant 4D energy equal to
$20$.

In Fig. 1, we illustrate the new 4D-band-limited ESS (4D-BL-ESS) scheme
through a trellis using 64-QAM that has 4 amplitudes in 1D:
$\left\\{1,3,5,7\right\\}$, a block length $N=8$, and a maximum accumulated
energy limit $E_{\mathrm{max}}=48$. $n$ and $L$ are respectively the column
and row indices of the nodes in the trellis. Each path in the trellis
correspond to a possible amplitude sequence. In each node, the number in black
$e$ is the intermediate accumulated energy up to that state, and the number in
red $T^{e}_{n}$ is the number of paths that lead from this state to a final
state at $n=N$. The accumulated 1D energy in the $4$-th and $8$-th positions
(red nodes) corresponds to the accumulated 4D energy in the $1$-st and $2$-nd
time slots respectively. Sequences in grey are deleted based on a 4D-energy
constraint instead of a 1D-constraint as initially done in [5]. For this
example, the constraint is: $4+28(i-1)\leq E_{4D,i}\leq 28+28(i-1)$ where
$E_{4D,i}$ is the accumulated 4D energy in $i$-th time slot. The construction
used in the toy example can be generalized by tuning the constraints on the
accumulated energy:

$A+iK_{1}\leq E_{4D,i}\leq\mathrm{min}(B+iK_{2},E_{\mathrm{max}})$ (1)

where the parameters $A,B,K_{1},K_{2}$ define the upper and lower limits of 4D
energies in the trellis.

Figure 2: Left: Electrical SNR versus Power per channel in the non-linear
channel. Right: GMI versus Power per channel. Simulation of conventional ESS,
K-ESS, 1D-BL-ESS, 4D-BL-ESS(Linear), and 4D-BL-ESS(Nonlinear) in $50$ GBaud
over $1\times 205$ km SSMF link, 5 WDM channels, 55 GHz channel spacing, net
rate per channel $400$ Gbit/s.

For ease of comparison between different schemes, we evaluate several metrics
related to the linear and non-linear performance using the same signaling
parameters as in [4, 5], i.e. a shaping rate of $1.5$ bits/amplitude, a block
length $N=108$, an LDPC FEC code with a rate of $r_{c}=5/6$ and acheive a net
bit rate of $8$ bits/4D symbol which is calculated as in [8]. For conventional
ESS, the net bit rate for each amplitude is determined only by the chosen
value of $E_{\mathrm{max}}$. However, for all three band-limited ESS and the
K-ESS, we start with a larger energy limit $E_{\mathrm{max}}$, hence with a
larger bit rate, and we then reduce the bit rate through constraints such as
band limiting the enrergy or limiting the kurtosis until we reach the same net
rate as the one of conventional ESS. In [4, 5], the average energy of 1D
amplitudes $\overline{E}_{1D}$ is used to assess the linear penalty, also the
variance of 1D amplitudes $\mathrm{\sigma^{2}_{1D}}$ and the kurtosis of 2D
amplitudes $\mathrm{\mu_{4,2D}}$ are the used metrics to assess the non-linear
gain. Hence, in Table I, we compute these metrics for five schemes: the
conventional ESS [1], the kurtosis-ESS (K-ESS) [4], the 1D-BL-ESS [5] and two
proposed 4D-BL schemes. We also add two 4D metrics: the energy dispersion
index (EDI) defined in [9], and the windowed kurtosis defined in [10]. These
two 4D metrics were shown to be better suited for assessing non-linear gains
of non-iid (independent and identically distributed) signaling schemes such as
short block-length PAS modulations [10, 9, 11]. The length of the averaging
windows $W_{1}$ for the EDI and $W_{2}$ for the windowed kurtosis are computed
as Eq. (2) in [9] and Eq. (7) in [10]. Their values in Table I are computed
for the first transmission scenario over a single-span $205$-km SSMF link
simulated in the next section. In this paper, we design and study two 4D-BL
schemes: the first one is optimized for the linear regime, is denoted 4D-BL-
ESS (Linear) and has a lower accumulated energy threshold $E_{\mathrm{max}}$
than its 1D-BL counterpart, thus a lower rate loss (linear gain) and almost
the same energy variations as 1D-BL-ESS; and the second scheme is optimized
for the non-linear regime, is denoted 4D-BL-ESS (Nonlinear) and has the same
accumulated energy threshold $E_{\mathrm{max}}$ as 1D-BL-ESS (hence, the same
rate loss too) but it has lower 4D energy variations. The rate loss for each
scheme is computed and reported in the last row of Table I. We can see that
all band-limited ESS schemes show lower energy variations than the
conventional ESS (first column), but have a higher rate loss due to the
addition of energy constraints that reduce the number of allowed amplitude
combinations. At the same net rate, K-ESS has the lowest 1D energy variation,
but the highest 4D energy variation, which makes its nonlinear performance
worse than the other schemes as will be shown in the next section.

Figure 3: Simulation of conventional ESS, 1D-BL-ESS, 4D-BL-ESS(Linear), and
4D-BL-ESS(Nonlinear) for a $110$ GBaud transmission over $5\times 80$ km SSMF
link, 5 WDM channels, $112.5$ GHz channel spacing, net rate per channel: $838$
Gbit/s. Electrical SNR (left) and AIR (middle) vs block length at optimal
power ($3$ dBm). EDI (right) with $W_{1}=200$ versus block length.

##  III Simulation setup and results

###  III-A Single-span $205$-km $50$GBaud transmission

We evaluate the performance of the two proposed 4D-BL-ESS with $N=108$ under
the same conditions as the ones used in [4, 5]. The transmitted signals are
$5$ dual-polarization WDM channels, with a baud rate of $50$ GBaud per
channel. The raw data rate is $600$ Gbit/s. For all ESS schemes, we use an
LDPC code of length $n=64800$ bits from the DVB-S2 standard with a rate of
$r_{c}=5/6$. The achieved net bit rate is $8$ bits per 4D symbol. The channel
spacing is $55$ GHz. Root-raised-cosine pulse shaping is applied with a roll-
off factor of $0.1$. The signal propagation is simulated through the split-
step Fourier method based on the nonlinear Schrödinger and Manakov equations.
SSMF fiber segments are simulated with a chromatic dispersion coefficient
$D=17\mathrm{ps/nm/km}$, non-linearity coefficient $\gamma=1.3\mathrm{(W\cdot
km)^{-1}}$, and an attenuation $\alpha=0.2~{}\mathrm{dB/km}$ at
$\lambda=1550~{}\mathrm{nm}$. Polarization mode dispersion (PMD) is added with
$0.04~{}\mathrm{ps/\sqrt{km}}$. The link is a single span of
$205~{}\mathrm{km}$. No laser phase noise is added at this point. After
propagation, the central channel is filtered. The signal undergoes matched
filtering, chromatic dispersion compensation, and genie-aided multiple-input-
multiple-output (MIMO) channel equalization. Finally, we use a fully-data-
aided phase filter with an averaging window of $64$ symbols to compensate for
the phase rotation induced by cross-phase modulation (XPM), as described in
[12] and we compute the electric SNR and the achievable information rate (AIR)
from the equalized signal. The electric signal-to-noise ratio
$\mathrm{SNR_{elec}}$ is the SNR computed from the constellation at the end of
the DSP chain at the receiver side (i.e., at the input of the decision
circuit) and the AIR is computed as Eq. (10) in [13]. For each launch power,
five different channel realizations are simulated.

In the left part of Fig. 2, we show the electric SNR versus the launch power
per channel for the five ESS schemes. We see that the 1D-BL, 4D-BL (Linear),
and 4D-BL (Nonlinear) ESS show respectively a $1$ dB, $0.9$ dB, and $1.1$ dB
increase in SNR at a power of $10$ dBm compared to conventional ESS, and a $1$
dB enhancement in optimal launch power. In the right part of Fig. 2, we show
the results in terms of AIR computed as Eq. (1) in [14] per 4D symbol versus
the optical power. The 4D-BL (linear) and (Nonlinear) show the best AIR. At a
launch power of $10$ dBm, the bit rate is increased by $0.07$ bit/4D compared
to the 1D-BL-ESS, and by $0.19$ bit/4D compared to conventional ESS. 4D-BL-ESS
(Nonlinear) also shows $0.13$ bit/4D gain compared to the 1D-BL-ESS in highly
non-linear regime. While K-ESS shows the lowest 1D energy variation as can be
seen from Table I, its 4D energy variation is higher than those of the BL-ESS
schemes, which is compliant with the measured SNR and AIR performance in Fig.
2.

###  III-B $5\times 80$-km multi-span $110$GBaud transmission

Next, we compare four schemes: conventional, 1D-BL and 4D-BL ESS (Linear and
Nonlinear) in a new scenario at a higher baud rate and over a multi-span link
configuration for different shaping block lengths $N$, to assess the scaling
of the linear and non-linear gains. We consider a $110$ GBaud transmission
over $5\times 80$ km SSMF link, $5$ WDM channels with a $112.5$ GHz channel
spacing. We also consider a laser linewidth of $100$ kHz at each of the
transmitter side and receiver side. The net rate per channel is $838$ Gbit/s
($880$ Gbit/s with $5\%$ pilot overhead). We compensate the added phase noise
using a pilot-aided phase estimation [15] with a $5\%$ overhead. The other DSP
blocks are the same as the ones used in the previous section. We study the
impact of the block length $N$ on the performance in terms of electric SNR and
AIR, both measured at the optimal launch power of $3$ dBm.

In Fig. 3, we show the variation of the electric SNR, the AIR and the EDI with
respect to the shaping block length $N=\\{60,108,200,300\\}$. For each block
length, five different channel realizations are simulated. First, we can seen,
by comparing the leftmost and the rightmost figures, that the measured SNR is
inversely correlated to EDI changes. However, we also notice that EDI only
gives a rough prediction of the performance as it fails to capture small
differences such as the gap between the performance of 1D-BL-ESS and 4D-BL-ESS
(Linear). Second, we can see that, at a short block length $N=60$, all BL-ESS
schemes fail to achieve a significant gain in SNR as the non-linear gains from
energy limitation are still small compared to the conventional ESS. Third,
when the rate loss decreases as the block length increases, and consequently
the linear penalty decreases, the energy fluctuations of conventional ESS
reduce the SNR while the SNR of the BL-ESS schemes remain stable up to
$N=300$. We did not increase the block length beyond $300$ to maintain a low
complexity for the distribution matching. Finally, if we assess the AIR gains
through the middle figure, we see that a net performance gain of BL-ESS
schemes with respect to conventional ESS is observed beyond a block length of
$200$. However, all the AIR values achieved by the BL-ESS schemes remain lower
than the AIR of conventional ESS for $N\leq 200$ due to their higher rate
loss. Our explanation of the limited gains of the BL-ESS for longer distances
and higher baud rates is the following: despite the optimization of the signal
characteristics at the transmitter side for low NLI distortions through the
addition of energy constraints, the chromatic dispersion gradually destroys
these characteristics during the transmission. Hence, the band-limited ESS can
only obtain limited nonlinear gains over long-distance and high baud-rate
transmission (gains achieved over the first few spans), and the gain in AIR is
also limited due to the higher rate loss of the band-limited ESS.

##  IV Conclusion

We proposed a new design of ESS schemes that constrains energy variations in
higher dimensions (four dimensions in this work) and that can be tuned to
boost either linear or nonlinear gains. The new scheme achieves higher gains
than state-of-the-art ESS schemes and is particularly interesting for short-
distance links as the trade-off between the achieved linear and non-linear
gains becomes harder to meet for longer distances because of the accumulation
of chromatic dispersion that increases again the energy variations of the
optical signal. Besides, our simulation results show that windowed 4D metrics
of energy variations predict the non-linear performance of ESS schemes better
than 1D and 2D metrics; however, their prediction is not yet accurate enough
and the search for better 4D metrics is still an interesting problem to solve.

##  V Acknowledgements

This work was supported by Huawei Technologies France.

## References

  * [1] Y. C. Gültekin, W. J. van Houtum, A. G. Koppelaar, F. M. Willems, J. Wim _et al._ , “Enumerative sphere shaping for wireless communications with short packets,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 2, pp. 1098–1112, 2019.
  * [2] P. Schulte and G. Böcherer, “Constant composition distribution matching,” _IEEE Transactions on Information Theory_ , vol. 62, no. 1, pp. 430–434, 2015.
  * [3] T. Fehenberger, D. S. Millar, T. Koike-Akino, K. Kojima, K. Parsons, and H. Griesser, “Analysis of nonlinear fiber interactions for finite-length constant-composition sequences,” _Journal of Lightwave Technology_ , vol. 38, no. 2, pp. 457–465, 2019.
  * [4] Y. C. Gültekin, A. Alvarado, O. Vassilieva, I. Kim, P. Palacharla, C. M. Okonkwo, and F. M. Willems, “Kurtosis-limited sphere shaping for nonlinear interference noise reduction in optical channels,” _Journal of Lightwave Technology_ , vol. 40, no. 1, pp. 101–112, 2021.
  * [5] Gültekin, Yunus Can and Alvarado, Alex and Vassilieva, Olga and Kim, Inwoong and Palacharla, Paparao and Okonkwo, Chigo M and Willems, Frans MJ, “Mitigating nonlinear interference by limiting energy variations in sphere shaping,” in _Optical Fiber Communication Conference_. Optica Publishing Group, 2022, pp. Th3F–2.
  * [6] K. Wu, G. Liga, A. Sheikh, Y. C. Gültekin, F. M. Willems, and A. Alvarado, “List-encoding CCDM: A nonlinearity-tolerant shaper aided by energy dispersion index,” _Journal of Lightwave Technology_ , vol. 40, no. 4, pp. 1064–1071, 2022.
  * [7] P. Skvortcov, I. Phillips, W. Forysiak, T. Koike-Akino, K. Kojima, K. Parsons, and D. S. Millar, “Huffman-coded sphere shaping for extended-reach single-span links,” _IEEE Journal of Selected Topics in Quantum Electronics_ , vol. 27, no. 3, pp. 1–15, 2021.
  * [8] L. Schmalen, “Probabilistic constellation shaping: Challenges and opportunities for forward error correction,” in _2018 Optical Fiber Communications Conference and Exposition (OFC)_ , 2018, pp. 1–3.
  * [9] K. Wu, G. Liga, A. Sheikh, F. M. Willems, and A. Alvarado, “Temporal energy analysis of symbol sequences for fiber nonlinear interference modelling via energy dispersion index,” _Journal of Lightwave Technology_ , vol. 39, no. 18, pp. 5766–5782, 2021.
  * [10] J. Cho, X. Chen, G. Raybon, D. Che, E. Burrows, S. Olsson, and R. Tkach, “Shaping lightwaves in time and frequency for optical fiber communication,” _Nature communications_ , vol. 13, no. 1, p. 785, 2022.
  * [11] Y. C. Gültekin, K. Wu, and A. Alvarado, “Temporal properties of enumerative shaping: Autocorrelation and energy dispersion index,” in _Signal Processing in Photonic Communications_. Optica Publishing Group, 2022, pp. SpM2I–3.
  * [12] T. Fehenberger, M. P. Yankov, L. Barletta, and N. Hanik, “Compensation of xpm interference by blind tracking of the nonlinear phase in wdm systems with qam input,” in _2015 European Conference on Optical Communication (ECOC)_. IEEE, 2015, pp. 1–3.
  * [13] G. Böcherer, “Probabilistic signal shaping for bit-metric decoding,” in _2014 IEEE International Symposium on Information Theory_ , 2014, pp. 431–435.
  * [14] A. Amari, S. Goossens, Y. C. Gültekin, O. Vassilieva, I. Kim, T. Ikeuchi, C. Okonkwo, F. M. Willems, and A. Alvarado, “Enumerative sphere shaping for rate adaptation and reach increase in WDM transmission systems,” in _45th European Conference on Optical Communication (ECOC 2019)_. IET, 2019, pp. 1–4.
  * [15] M. S. Neves, A. Carena, A. Nespola, P. P. Monteiro, and F. P. Guiomar, “Joint carrier-phase estimation for digital subcarrier multiplexing systems with symbol-rate optimization,” _Journal of Lightwave Technology_ , vol. 39, no. 20, pp. 6403–6412, 2021.


# On the optimality of MCMC integration for functions with unbounded second
moment

Julian Hofstadler Faculty of Computer Science and Mathematics, University of
Passau, Innstraße 33, 94032 Passau, Germany.
Email<EMAIL_ADDRESS>
###### Abstract

We study the Markov chain Monte Carlo (MCMC) estimator for numerical
integration for functions that do not need to be square integrable w.r.t. the
invariant distribution. For chains with a spectral gap we show that the
absolute mean error for $L^{p}$ functions, with $p\in(1,2)$, decreases like
$n^{1/p-1}$, which is known to be the optimal rate. This improves currently
known results where an additional parameter $\delta>0$ appears and the
convergence is of order $n^{(1+\delta)/p-1}$.

Keywords: Markov chain Monte Carlo, spectral gap, absolute mean error

Classification: 65C05; 60J22; 65C20

##  1 Introduction

Let $(\mathcal{X},\mathcal{F}_{\mathcal{X}},\pi)$ be a probability space and
$f\colon\mathcal{X}\to\mathbb{R}$ measurable as well as $\pi$-integrable. For
a random variable $X\sim\pi$ we are interested in approximating the
expectation

$\mathbb{E}[f(X)]=\int_{\mathcal{X}}f(x)\pi({\rm d}x)=\pi(f).$

A common approach is to use an MCMC method. Requiring the density of $\pi$
only in non-normalised form, many MCMC algorithms provide powerful tools in
scientific and statistical applications. The main idea behind these approaches
is to construct a Markov chain $(X_{n})_{n\in\mathbb{N}}$, having $\pi$ as the
invariant distribution, and to estimate $\pi(f)$ via

$S_{n}f=\frac{1}{n}\sum_{j=1}^{n}f(X_{j}).$

Under fairly mild conditions one has $S_{n}f\to\pi(f)$ almost surely as
$n\to\infty$, cf. [AG11] or [MT12, Chapter 17]. This ensures (strong)
consistency of the MCMC estimator, yet, it is clearly of interest to have non-
asymptotic error bounds. For instance, given some $p\in[1,\infty)$, one may
consider the $p$-mean error

$\mathbb{E}\left[\left|S_{n}f-\pi(f)\right|^{p}\right],$ (1)

however, also other criteria are feasible, see e.g. [KR19].

Setting $p=2$ in (1) we speak of the mean squared error, and in different
settings explicit bounds are known, e.g. under a Wasserstein contraction
assumption [JO10], spectral gap conditions [Rud12], or if one has
geometric/polynomial ergodicity [ŁMN13]. For the mean squared error to make
sense we require a finite second moment of $f$, i.e.
$\left\|f\right\|_{L^{2}(\pi)}^{2}=\pi(f^{2})<\infty$. On the other hand, the
absolute mean error, given by $\mathbb{E}|S_{n}f-\pi(f)|$, is well defined and
finite as long as $f$ is $\pi$-integrable. Bounds for the absolute mean error
for functions with $\left\|f\right\|_{L^{p}(\pi)}^{p}=\pi(|f|^{p})<\infty$,
where $p<2$, are still rare. In [RS15] it is shown that under a spectral gap
condition for any $p\in(1,2)$ one has

$\sup_{\left\|f\right\|_{L^{p}(\pi)}\leq
1}\mathbb{E}[|S_{n}f-\pi(f)|]\leq\frac{C}{n^{1-(1+\delta)/p}},$ (2)

with constants $\delta>0$ and $C\in(0,\infty)$.

Proposition 1 of Section 2.2.9 (with $k=0$) in [Nov88], see also [Hei94,
Section 5], shows that in general we have the following lower bound

$\sup_{\left\|f\right\|_{L^{p}(\pi)}\leq
1}\mathbb{E}[|S_{n}f-\pi(f)|]\geq\frac{c}{n^{1-1/p}},$ (3)

where $c\in(0,\infty)$ is a constant independent of $n$.

Even though $\delta>0$ in (2) may be chosen arbitrarily small, it is natural
to ask whether it can be removed completely, such that one would have the same
rate as in (3). Under the (strong) assumption of uniform ergodicity and
reversibility we know that this is the case, cf. [RS15, Theorem 1]. To the
best of our knowledge this is the only situation where optimal rates are
known. The goal of this note is to extend that result to the spectral gap
setting, see Corollary 2.4, where we show

$\sup_{\left\|f\right\|_{L^{p}(\pi)}\leq
1}\mathbb{E}[|S_{n}f-\pi(f)|]\leq\frac{\widetilde{C}}{n^{1-1/p}},$

with an explicit expression for the constant $\widetilde{C}$.

Let us sketch the proof: The main idea is to employ the Riesz-Thorin
interpolation theorem, a technique which goes back at least to [Hei94] in
Monte Carlo theory and was also used to derive (2) in [RS15]. In contrast to
[RS15] we first show a bound for $L^{2}$ functions, see Proposition 3.2, which
leads to the optimal rate after interpolation. It is worth mentioning that
Proposition 3.2 extends results of [Rud12], see Remark 3.3.

The rest of this note is organised as follows: In Section 2 we state and
discuss our assumptions as well as main results. The proofs, together with
some intermediate results, can be found in Section 3.

##  2 Error bounds for MCMC integration

This section contains our main results, see Theorem 2.3 and Corollary 2.4,
together with all required notation and assumptions. Let us start with
specifying the general setting.

We assume that the state space $(\mathcal{X},\mathcal{F}_{\mathcal{X}})$ is
Polish with $\mathcal{F}_{\mathcal{X}}$ being countably generated. Let $K$ be
a Markov kernel and $\nu$ a probability measure, called the initial
distribution, both defined on $(\mathcal{X},\mathcal{F}_{\mathcal{X}})$. Then,
the Markov chain corresponding to $K$ and $\nu$, say
$(X_{n})_{n\in\mathbb{N}}$, is defined on a probability space
$(\Omega,\mathcal{F},\mathbb{P}_{\nu})$. (In particular, such a probability
space exists.) We assume that $\pi$ is the unique invariant distribution of
$K$ and that $(X_{n})_{n\in\mathbb{N}}$ is $\psi$-irreducible. For definitions
and further details we refer to [MT12, DMPS18].

Given $p\in[1,\infty)$, we define $L^{p}(\pi)$ as the set of all (equivalence
classes111Two functions are considered as equivalent if they are equal $\pi$
almost everywhere. of) functions $f\colon\mathcal{X}\to\mathbb{R}$ such that
$\left\|f\right\|_{L^{p}(\pi)}<\infty$. Similar, $L^{\infty}(\pi)$ denotes the
set of (equivalence classes of) functions with finite
$\left\|f\right\|_{L^{\infty}(\pi)}=\operatornamewithlimits{ess\,sup}_{x\in\mathcal{X}}|f(x)|$.
Recall that $L^{2}(\pi)$ (equipped with $\langle
f,g\rangle_{L^{2}(\pi)}=\int_{\mathcal{X}}f(x)g(x)\pi({\rm d}x)$) is a Hilbert
space and so is the following (closed) subspace

$L^{2}_{0}(\pi):=\left\\{f\in L^{2}(\pi)\colon\pi(f)=0\right\\}.$

Let $p\in[1,\infty]$, then the Markov kernel $K$ induces a linear operator
$K\colon L^{p}(\pi)\to L^{p}(\pi)$ via $f\mapsto
Kf(\cdot)=\int_{\mathcal{X}}f(x^{\prime})K(\cdot,{\rm d}x^{\prime})$. Indeed,
(the operator) $K$ is well defined and one has
$\left\|K\right\|_{L^{p}(\pi)\to L^{p}(\pi)}=1$, we refer to [Rud12, Section
3.1] for further details.

By $\operatorname{\text{Id}}$ we denote the identity on $L^{2}(\pi)$. The
following condition about the operator $\operatorname{\text{Id}}-K$,
restricted to $L^{2}_{0}(\pi)$, is our main assumption.

######  Assumption 2.1.

Assume that $\operatorname{\text{Id}}-K$, considered as operator from
$L^{2}_{0}(\pi)$ to $L^{2}_{0}(\pi)$, has a linear and bounded inverse with

$\left\|(\operatorname{\text{Id}}-K)^{-1}\right\|_{L^{2}_{0}(\pi)\to
L^{2}_{0}(\pi)}\leq s<\infty.$

######  Remark 2.2.

The invertibility of $\operatorname{\text{Id}}-K$, restricted to a suitable
subspace of $L^{2}(\pi)$, was also studied in [Mat99] for uniformly ergodic
chains and [Mat04] for $V$-uniformly ergodic chains. In particular, the
existence of $(\operatorname{\text{Id}}-K)^{-1}$ on an appropriate subspace
was used there to characterise the convergence behaviour of the mean squared
error.

We note that Assumption 2.1 is closely related to a spectral gap, there are,
however, different definitions for a spectral gap:

  * •

Some authors, see e.g. [DMPS18, ALPW23], say $K$ admits a(n) (absolute
$L^{2}$) spectral gap if $\sup_{\lambda\in\mathcal{S}_{0}}|\lambda|<1$, where
$\mathcal{S}_{0}$ is the spectrum of $K\colon L^{2}_{0}(\pi)\to
L^{2}_{0}(\pi)$. This is equivalent to the existence of some $m\in\mathbb{N}$
such that $\left\|K^{m}\right\|_{L^{2}_{0}(\pi)\to L^{2}_{0}(\pi)}<1$, cf.
[DMPS18, Proposition 22.2.4].

  * •

On the other hand, a different definition, for instance used in [Rud12, HSV14,
RS15, NRS21], is to say that $K$ admits a(n) (absolute $L^{2}$) spectral gap
if $\left\|K\right\|_{L^{2}_{0}(\pi)\to L^{2}_{0}(\pi)}<1$.

If $K$ is reversible (which implies that the corresponding Markov operator is
self-adjoint on $L^{2}(\pi)$), then both definitions are equivalent.

Let us emphasize that either of the above definitions of a spectral gap
implies that Assumption 2.1 is true, also in the non-reversible case. Spectral
gap results were established for a number of MCMC methods, see for instance
[HSV14, NRS21, ALPW23], see also [Rud12, Section 3.4] and [RR97, Theorem 2.1].
Moreover, we note that under Assumption 2.1 we cover the setting of Theorems 1
and 2 of [RS15].

######  Theorem 2.3.

Let Assumption 2.1 be true, $p\in(1,2)$ and assume that $\nu$ is absolutely
continuous w.r.t. $\pi$ with Radon-Nikodym derivative $\frac{{\rm d}\nu}{{\rm
d}\pi}\in L^{\infty}(\pi)$. Then, for any $f\in L^{p}(\pi)$ and any
$n\in\mathbb{N}$ we have that

$\mathbb{E}_{\nu}\left[\left|S_{n}f-\pi(f)\right|^{p}\right]\leq{C_{1}}\frac{\left\|f\right\|_{L^{p}(\pi)}^{p}}{n^{{p-1}}},$

where $C_{1}=\left(1+\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}\right)^{2-p}\left(12s^{2}\left\|\frac{{\rm
d}\nu}{{\rm d}\pi}\right\|_{L^{\infty}(\pi)}\right)^{p-1}$.

The former results together with Hoelder’s inequality implies the following
estimate for the absolute mean error where the rate of convergence is optimal.

######  Corollary 2.4.

Let Assumption 2.1 be true, $p\in(1,2)$ and assume that $\nu$ is absolutely
continuous w.r.t. $\pi$ with Radon-Nikodym derivative $\frac{{\rm d}\nu}{{\rm
d}\pi}\in L^{\infty}(\pi)$. Then, for any $f\in L^{p}(\pi)$ and any
$n\in\mathbb{N}$ we have that

$\mathbb{E}_{\nu}\left[\left|S_{n}f-\pi(f)\right|\right]\leq{C_{2}}\frac{\left\|f\right\|_{L^{p}(\pi)}}{n^{{1-\frac{1}{p}}}},$

where $C_{2}=\left(1+\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}\right)^{2/p-1}\left(12s^{2}\left\|\frac{{\rm
d}\nu}{{\rm d}\pi}\right\|_{L^{\infty}(\pi)}\right)^{1-1/p}$.

######  Remark 2.5.

We note that Theorem 2.3 and Corollary 2.4 are still true (with the same
bounds) if we replace $S_{n}f$ by
$S_{n,n_{0}}f=\frac{1}{n}\sum_{j=1}^{n}f(X_{j+n_{0}})$, which corresponds to
using a burn in of length $n_{0}\in\mathbb{N}$.

######  Remark 2.6.

In Theorem 2.3 and Corollary 2.4 the quantity $\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}$ appears. In some sense this penalises our
choice of the initial distribution $\nu$ which is allowed to differ from the
target $\pi$. In the setting where a fixed computational budget is available
it may be worth spending some effort to find a “good” initial distribution
$\nu$. However, discussing optimal choices of $\nu$ w.r.t. different
theoretical and/or practical aspects is beyond the scope of this note.

##  3 Proofs

In this section we prove our main results. We start with an auxiliary lemma
which relates the Monte Carlo sum $S_{n}f$ to a martingale. To do so we use
the invertibility of $\operatorname{\text{Id}}-K$ to compute solutions of
Poisson’s equation. See also [MT12, Chapter 17.4.1] for further details about
that topic.

For any $n\in\mathbb{N}$ we set $\mathcal{F}_{n}=\sigma(X_{j}\colon 1\leq
j\leq n+1)$.

######  Lemma 3.1.

Let Assumption 2.1 hold true and let $h\in L^{2}_{0}(\pi)$. Then, setting
$\hat{h}=(\operatorname{\text{Id}}-K)^{-1}h$, we have that
$\mathbb{P}_{\nu}$-almost surely

$S_{n}h=\frac{1}{n}M_{n}+\frac{1}{n}\left(\hat{h}(X_{1})-\hat{h}(X_{n+1})\right),$

with $M_{n}=\sum_{j=1}^{n}\left(\hat{h}(X_{j+1})-K\hat{h}(X_{j})\right)$ being
an $(\mathcal{F}_{n})_{n\in\mathbb{N}}$ martingale. For its differences
$\Delta_{j}=\hat{h}(X_{j+1})-K\hat{h}(X_{j})$ we have

$\forall
j\in\mathbb{N}\colon\quad\mathbb{E}_{\nu}\left[|\Delta_{j}|^{2}\right]\leq
2s^{2}\left\|h\right\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}.$

###### Proof.

Set $T=\operatorname{\text{Id}}-K$. Note that by assumption $\hat{h}\in
L^{2}_{0}(\pi)$ and

$h=TT^{-1}h=T\hat{h}=\hat{h}-K\hat{h}.$

Thus $\mathbb{P}_{\nu}$-almost surely

$\displaystyle S_{n}h$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\hat{h}(X_{j})-K\hat{h}(X_{j})$
$\displaystyle=\frac{1}{n}\sum_{j=1}^{n}\left(\hat{h}(X_{j+1})-K\hat{h}(X_{j})\right)+\frac{1}{n}\left(\hat{h}(X_{1})-\hat{h}(X_{n+1})\right)$
$\displaystyle=\frac{1}{n}M_{n}+\frac{1}{n}\left(\hat{h}(X_{1})-\hat{h}(X_{n+1})\right).$

Since $(X_{n})_{n\in\mathbb{N}}$ is a Markov chain with kernel $K$ for any
$n\geq 2$ holds

$\mathbb{E}_{\nu}[M_{n}|\mathcal{F}_{n-1}]=M_{n-1}+\mathbb{E}_{\nu}\left[\hat{h}(X_{n+1})-K\hat{h}(X_{n})\;\big{|}\;\mathcal{F}_{n-1}\right]=M_{n-1}.$

It remains to show the bound for the differences $\Delta_{j}$, which also
implies that $\mathbb{E}_{\nu}|M_{n}|<\infty$ for any $n\in\mathbb{N}$ and
hence that $(M_{n})_{n\in\mathbb{N}}$ is a martingale. Thereto observe that
$\|\hat{h}\|_{L^{2}(\pi)}\leq s\left\|h\right\|_{L^{2}(\pi)}$ and that
$\left\|K\right\|_{L^{1}(\pi)\to L^{1}(\pi)}=1$. Hence, Hoelder’s inequality
shows that for any $j\in\mathbb{N}$ holds

$\displaystyle\mathbb{E}_{\nu}\left[|\hat{h}(X_{j+1})|^{2}\right]=\int_{\mathcal{X}}K^{j}|\hat{h}|^{2}(x)\frac{{\rm
d}\nu}{{\rm d}\pi}(x)\pi({\rm
d}x)\leq\left\|K^{j}|\hat{h}|^{2}\right\|_{L^{1}(\pi)}\left\|\frac{{\rm
d}\nu}{{\rm d}\pi}\right\|_{L^{\infty}(\pi)}$
$\displaystyle\leq\|\hat{h}\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}\leq
s^{2}\|h\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}.$

Similarly,
$\mathbb{E}_{\nu}[(K\hat{h}(X_{j}))^{2}]\leq\mathbb{E}_{\nu}[K\hat{h}^{2}(X_{j})]\leq
s^{2}\|h\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}$. The bound for
$\mathbb{E}_{\nu}|\Delta_{j}|^{2}$ now follows by using the inequality
$(\alpha+\beta)^{2}\leq 2(\alpha^{2}+\beta^{2})$, being valid for any
$\alpha,\beta\in\mathbb{R}$. ∎

Based on the above martingale decomposition of Lemma 3.1 we now show a bound
for the mean squared error.

######  Proposition 3.2.

Let Assumption 2.1 be true and assume that $\nu$ is absolutely continuous
w.r.t. $\pi$ with Radon-Nikodym derivative $\frac{{\rm d}\nu}{{\rm d}\pi}\in
L^{\infty}(\pi)$. Then, for any $f\in L^{2}(\pi)$ we have that

$\mathbb{E}_{\nu}\left[\left|S_{n}f-\pi(f)\right|^{2}\right]\leq
6s^{2}\left(1+\frac{1}{n}\right)\frac{\left\|f\right\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm
d}\nu}{{\rm d}\pi}\right\|_{L^{\infty}(\pi)}}{n}.$

###### Proof.

Let $h=f-\pi(f)\in L^{2}_{0}(\pi)$. Lemma 3.1 and the inequality
$(\alpha+\beta+\gamma)^{2}\leq 3(\alpha^{2}+\beta^{2}+\gamma^{2})$, which is
true for any real numbers $\alpha,\beta,\gamma$, yield

$\mathbb{E}_{\nu}\left[|S_{n}f-\pi(f)|^{2}\right]=\mathbb{E}_{\nu}\left[(S_{n}h)^{2}\right]\leq\frac{3}{n^{2}}\mathbb{E}_{\nu}\left[M_{n}^{2}\right]+\frac{3}{n^{2}}\mathbb{E}_{\nu}[\hat{h}(X_{1})^{2}+\hat{h}(X_{n+1})^{2}],$

with $M_{n},\hat{h}$ as specified in Lemma 3.1. By the martingale property
(and again using the notation/results of Lemma 3.1) it follows that

$\mathbb{E}_{\nu}\left[M_{n}^{2}\right]=\mathbb{E}_{\nu}\left[\sum_{j=1}^{n}\Delta_{j}^{2}\right]\leq\sum_{j=1}^{n}\mathbb{E}_{\nu}\left[\Delta_{j}^{2}\right]\leq
2ns^{2}\|h\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}.$

Furthermore, the same arguments as in the proof of Lemma 3.1 reveal

$\mathbb{E}_{\nu}[\hat{h}(X_{1})^{2}+\hat{h}(X_{n+1})^{2}]\leq
2s^{2}\|h\|_{L^{2}(\pi)}^{2}\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}.$

Combining these estimates with
$\left\|h\right\|_{L^{2}(\pi)}=\left\|f\right\|_{L^{2}(\pi)}^{2}-\pi(f)^{2}\leq\left\|f\right\|_{L^{2}(\pi)}^{2}$
yield the desired result. ∎

######  Remark 3.3.

Note that Proposition 3.2 can be seen as extension of two results in [Rud12]:

  * •

It extends Theorem 3.34 to the case where $K$ is not required to be
reversible222Note that $L^{1}$-exponential convergence implies
$L^{2}$-exponential convergence, e.g. by an application of the Riesz-Thorin
interpolation theorem, see [Gra14, Theorem 1.3.4]. Thus Assumption 2.1 holds
true in this case..

  * •

It extends Theorem 3.41 to the $L^{2}$ setting.

Finally, we turn to the proof of our main result.

###### Proof of Theorem 2.3.

We write $L^{q}(\mathbb{P}_{\nu})$ for the set of random variables $X$ with
finite
$\left\|X\right\|_{L^{q}(\mathbb{P}_{\nu})}^{q}=\mathbb{E}_{\nu}|X|^{q}$. For
$n\in\mathbb{N}$ consider the linear mapping $T_{n}f=S_{n}f-\pi(f)$ and note
that

$\mathbb{E}_{\nu}|T_{n}f|\leq\frac{1}{n}\sum_{j=1}^{n}\mathbb{E}_{\nu}|f(X_{j})|+\pi(|f|)=\frac{1}{n}\sum_{j=1}^{n}\int_{\mathcal{X}}K^{j-1}|f|(x)\frac{{\rm
d}\nu}{{\rm d}\pi}(x)\pi({\rm d}x)+\pi(|f|).$

By Hoelder’s inequality and the fact that $\left\|K\right\|_{L^{1}(\pi)\to
L^{1}(\pi)}=1$ it follows that $\left\|T_{n}\right\|_{L^{1}(\pi)\to
L^{1}(\mathbb{P}_{\nu})}\leq 1+\left\|\frac{{\rm d}\nu}{{\rm
d}\pi}\right\|_{L^{\infty}(\pi)}$.

Additionally, by Proposition 3.2 holds

$\left\|T_{n}\right\|_{L^{2}(\pi)\to
L^{2}(\mathbb{P}_{\nu})}\leq\frac{1}{\sqrt{n}}\left(12s^{2}\left\|\frac{{\rm
d}\nu}{{\rm d}\pi}\right\|_{L^{\infty}(\pi)}\right)^{1/2}.$

An application of the Riesz-Thorin interpolation theorem, see [Gra14, Theorem
1.3.4] in the setting $p_{0}=q_{0}=1$; $p_{1}=q_{1}=2$ and $\theta=2-2/p$,
yields that

$\left\|T_{n}\right\|_{L^{p}(\pi)\to
L^{p}(\mathbb{P}_{\nu})}\leq\left\|T_{n}\right\|_{L^{1}(\pi)\to
L^{1}(\mathbb{P}_{\nu})}^{1-\theta}\left\|T_{n}\right\|_{L^{2}(\pi)\to
L^{2}(\mathbb{P}_{\nu})}^{\theta},$

which, together with the above estimates, implies the result. ∎

## Acknowledgements

The author would like to thank Mareike Hasenpflug and Daniel Rudolf for
valuable discussions which helped to significantly improve the presentation of
the paper. JH is supported by the DFG within project 432680300 – SFB 1456
subproject B02.

## References

  * [AG11] Søren Asmussen and Peter W. Glynn, _A new proof of convergence of MCMC via the ergodic theorem_ , Statistics & Probability Letters 81 (2011), no. 10, 1482–1485.
  * [ALPW23] Christophe Andrieu, Anthony Lee, Sam Power, and Andi Q. Wang, _Explicit convergence bounds for Metropolis Markov chains: isoperimetry, spectral gaps and profiles_ , arXiv preprint arXiv:2211.08959, 2023.
  * [DMPS18] Randal Douc, Éric Moulines, Pierre Priouret, and Philippe Soulier, _Markov chains_ , Springer, 2018.
  * [Gra14] Loukas Grafakos, _Classical Fourier analysis_ , third ed., Graduate Texts in Mathematics, vol. 249, Springer, 2014.
  * [Hei94] Stefan Heinrich, _Random approximation in numerical analysis_ , Proceedings of the Conference “Functional Analysis” Essen (1994), 123–171.
  * [HSV14] Martin Hairer, Andrew M. Stuart, and Sebastian J. Vollmer, _Spectral gaps for a Metropolis–Hastings algorithm in infinite dimensions_ , The Annals of Applied Probability 24 (2014), no. 6, 2455–2490.
  * [JO10] Aldéric Joulin and Yann Ollivier, _Curvature, concentration and error estimates for Markov chain Monte Carlo_ , The Annals of Probability 38 (2010), no. 6, 2418 – 2442.
  * [KR19] Robert J. Kunsch and Daniel Rudolf, _Optimal confidence for Monte Carlo integration of smooth functions_ , Advances in Computational Mathematics 45 (2019), no. 5, 3095–3122.
  * [ŁMN13] Krzysztof Łatuszyński, Błażej Miasojedow, and Wojciech Niemiro, _Nonasymptotic bounds on the estimation error of MCMC algorithms_ , Bernoulli 19 (2013), no. 5A, 2033–2066.
  * [Mat99] Peter Mathé, _Numerical integration using Markov Chains_ , Monte Carlo Methods and Applications 5 (1999), no. 4, 325–343.
  * [Mat04] , _Numerical integration using v-uniformly ergodic Markov Chains_ , Journal of Applied Probability 41 (2004), no. 4, 1104–1112.
  * [MT12] Sean P. Meyn and Richard L. Tweedie, _Markov chains and stochastic stability_ , second ed., Cambridge University Press, 2012.
  * [Nov88] Erich Novak, _Deterministic and stochastic error bounds in numerical analysis_ , Lecture notes in Mathematics, vol. 1349, Springer, 1988.
  * [NRS21] Viacheslav Natarovskii, Daniel Rudolf, and Björn Sprungk, _Quantitative spectral gap estimate and Wasserstein contraction of simple slice sampling_ , The Annals of Applied Probability 31 (2021), no. 2, 806 – 825.
  * [RR97] Gareth O. Roberts and Jeffrey S. Rosenthal, _Geometric ergodicity and hybrid Markov chains_ , Electronic Communications in Probability 2 (1997), no. none, 13 – 25.
  * [RS15] Daniel Rudolf and Nikolaus Schweizer, _Error bounds of MCMC for functions with unbounded stationary variance_ , Statistics & Probability Letters 99 (2015), 6–12.
  * [Rud12] Daniel Rudolf, _Explicit error bounds for Markov chain Monte Carlo_ , Dissertationes Math. 485 (2012).


# Inferring Compensatory Kinase Networks in Yeast using Prolog

George A. Elder Queen Mary University of LondonResearch partially funded by
Medical Research Council MRCSchool of Biological and Behavioural Sciences
QMUL
London, United Kingdom<EMAIL_ADDRESS>Queen Mary University of
LondonSchool of Biological and Behavioural Sciences
QMUL
London, United Kingdom Conrad Bessant Queen Mary University of LondonSchool of
Biological and Behavioural Sciences
QMUL
London, United Kingdom<EMAIL_ADDRESS>
###### Abstract

Signalling pathways are conserved across different species, therefore making
yeast a model organism to study these via disruption of kinase activity. Yeast
has 159 genes that encode protein kinases and phosphatases, and 136 of these
have counterparts in humans. Therefore any insight in this model organism
could potentially offer indications of mechanisms of action in the human
kinome. This study utilises a Prolog-based approach, data from a yeast kinase
deletions strains study and publicly available kinase-protein associations.
Prolog, a programming language that is well-suited for symbolic reasoning is
used to reason over the data and infer compensatory kinase networks. This
approach is based on the idea that when a kinase is knocked out, other kinases
may compensate for this loss of activity. Background knowledge on kinases
targeting proteins is used to guide the analysis. This knowledge is used to
infer the potential compensatory interactions between kinases based on the
changes in phosphorylation observed in the phosphoproteomics data from the
yeast study. The results demonstrate the effectiveness of the Prolog-based
approach in analysing complex cell signalling mechanisms in yeast. The
inferred compensatory kinase networks provide new insights into the regulation
of cell signalling in yeast and may aid in the identification of potential
therapeutic targets for modulating signalling pathways in yeast and other
organisms.

##  1 Introduction

As with any organism, so with Saccharomyces cerevisiae, protein
phosphorylation plays a crucial role in various signalling pathways that
control cellular activities such as gene expression, protein synthesis, DNA
replication, RNA processing, cell cycle regulation, metabolism, transport of
vesicles and other organelles, response to environmental stress and nutrients,
and cell differentiation. The specific actions of kinases and substrates are
largely determined by protein interactions, which are effected via
phosphorylation [11].

It is also known that many of these signalling pathways are conserved across
different species, so yeast can be used as a model organism to study
phosphorylation networks. Yeast has 159 genes that encode protein kinases and
phosphatases, and 136 of these have counterparts in humans [9]. Therefore any
insight in this model organism could potentially offer indications of
mechanisms of action in the human kinome.

Ours is not the first application of symbolic AI in the study of yeast and its
pathways more specifically. A major inspiration for applying techniques that
fall under the symbolic framework came from the extensive work done by Prof
R.King. Particularly, work of the type found in this study by A.Clare and
R.King [3] showed that it is possible to build a knowledge base via the use of
a declarative logic programming language. The seminal work on Robot Scientist
Project also offered insight in how logic programming (in this case inductive)
can be used to model the metabolism of yeast and how deletion in associated
genes can affect their respective processes [7]. Further use of a Prolog based
approach can be found in [13] where a Flux Balance Analysis model was combined
with genome related information from the KEGG database to predict the growth
patterns of yeast. More recently, network reconstruction as well as new
structure hypotheses were made using ProbLog, with the aim of providing new
insights to experimental biologists [5].

However, all of the above approaches have focused on the use of gene-based
background knowledge and experimentally derived facts to reach their
conclusions. These lack the granularity of a proteomics and even more so a
phosphoproteomics based approach.

##  2 Methods

###  2.1 Toy Model development

Initially a Toy Model version of our logic program was developed, in order to
test the feasibility and applicability of our proposed methodology. The
perturbations modelled were gene deletions and the question thus became
whether silencing (or deletion) of a given kinase has actually worked. Within
the context of the study chosen (to be described in detail in subsequent
section), specific gene deletions were confirmed by either proteomics or PCR
assays. Therefore and due to the nature of gene deletions in yeast it was
unnecessary to identify off targets or known targets of a given perturbation.

In Figure 1 the overall structure of our yeast Toy Model can be seen. Edges
represent either perturbations (deletions in this case), observational facts,
presence of a phosphosite on a kinase and kinase phosphosite associations.
Nodes represent either phosphosites, perturbations, kinases or other proteins.

Figure 1: A Toy Model illustrates the minimal number of interactions needed to
understand a biological scenario involving kinases (k), proteins (p), and
their phosphosites (s) in response to an initial perturbation ($\delta$). In
our model, the perturbation refers to a gene encoding a specific kinase being
deleted. Our background information in this implementation has levels of
granularity, in this case the kinase targeting a phosphosite on a protein
version is depicted.

To determine whether a kinase has been successfully silenced, development of
the logic program began with the first iteration of rules. When queried, it is
satisfied in the case where a kinase K (1), for which a known phosphosite S
target on a protein P is found (2). Lastly, (3) checks whether a corresponding
mass spectrometry fact following deletion of K can be matched. If all these
conditions are met then the interpreter responds with true thus confirming
that the $\delta$gene perturbation encoding kinase K has been successful.

⬇

1doesSilofKwork(deletionOF((kinase(K))) :-

2 phosphorylates(kinase(K), phosphosite(S), protein(P)),

3 perturbs(deletionOF((kinase(K))),phosphosite(S), occupancy(down),_).

In the next rule, developed in an iterative manner, a “uniqueness check” was
added. This establishes whether a phosphosite S is uniquely targeted by a
kinase K. This is achieved by establishing a rule and checking for its
negation as part of this second iteration. In Prolog notation:

⬇

1doesSilofKwork2(deletionOF((kinase(K))) :-

2 phosphorylates(kinase(K), phosphosite(S), protein(P)),

3 \\+ sharedtarget(kinase(K), kinase(_K2), phosphosite(S)),

4 perturbs(deletionOF((kinase(K))),phosphosite(S), occupancy(down),_).

The rule sharedtarget(kinase(K), kinase(_K2), phosphosite(S)) includes two
goals matching phosphorylation facts on shared phosphosites then checking that
the two phosphorylating are not the same.

⬇

1sharedtarget(kinase(K), kinase(K2), phosphosite(S)) :-

2 (phosphorylates(kinase(K), phosphosite(S), kinase(P)),

3 phosphorylates(kinase(K2), phosphosite(S), kinase(P)),

4 K1 \= K2).

In the toy model, the responses to these two iterations were as expected.
However, there is still a concern that when applied to the larger dataset with
more than five perturbation-affected phosphosites, the true responses will be
based on only one phosphosite instance. Therefore the addition of a majority
check was decided. The final version of the rule took the following form:

⬇

1doesSilofKwork3(deletionOF((kinase(K))) :-

2 phosphorylates(kinase(K), phosphosite(S), kinase(P),

3 \\+ sharedtarget(kinase(K), kinase(_K2), phosphosite(S)),

4 findall((phosphosite(S),

5 perturbs(deletionOF((kinase(K))),phosphosite(S), occupancy(up),_),

6 perturbs(deletionOF((kinase(K))),phosphosite(S), occupancy(unaffected),_),

7 List1)),

8 findall((phosphosite(S),

9 perturbs(deletionOF((kinase(K))),phosphosite(S), occupancy(down),_),

10 List2)))),

11 length(List1, Length1), length(List2, Length2),

12 compare( ’>’, Length2, Length1).

This rule adds the majority check to the previous rule version. This creates
two lists, one containing all known targets of the kinase that are either “up”
or “unaffected” and another that contains all the ones described as “down”,
then compares the length of the two. The aim of this is to establish whether
the majority of the known targets were “down” rather than not, therefore
confirming the silencing of the kinase.

###  2.2 Experimental facts and background knowledge sources

A recently published study that looked at the effect of kinase perturbation at
the phosphoproteomic level was chosen as the source of our experimental data.
In [8], a systems-level proteomic and phosphoproteomic analysis was carried
out on 110 yeast single-kinase or phosphatase deletion strains under standard
growth conditions. The 110 were split in 84 kinases and 26 phosphatases. They
employed various methods, including traditional enrichment analysis,
$\delta$gene-$\delta$gene correlation networks, and molecular covariance
networks in order to analyse the functional relationships between these active
proteins.

In total their chosen deletions, which will be referred to as perturbations in
this section, cover 82% of all possible yeast and phosphatase deletion
strains. Through their experimental and mass spectrometry workflow they were
able to identify more than 4,600 and 13,000 proteins and phosphosites,
respectively. Of interest is the finding that they were able to capture a
large part of regulated phosphorylation events as well as 30% of those being
newly captured. This was entirely attributed to their ability to normalise
their data with protein abundance. This is also the data that we elected to
include in our ’perturbs/n’ predicate format.

Key constituents of our background knowledge base are the associations between
kinases and their targets both at the protein and phosphosite levels. For
this, two main sources were chosen, namely the Yeast Kinase Interaction
Database [12] (accessible at http://www.moseslab.csb.utoronto.ca/KID/) and the
Yeast Kinome database [2] (accessible at https://thebiogrid.org/project/2).

The Yeast Kinase Interaction Database (KID) is a resource that contains data
on phosphorylation events from various high- and low-throughput experiments.
It includes a total of literature-curated low-throughput (6,225) and high-
throughput (21,990) interactions, stemming from over 35,000 experiments. It
includes 517 high-confidence (or 853 low-confidence) kinase-substrate pairs
depending on the cutoffs used for the metric provided in the study. On the
other hand the Yeast Kinome database contains 1,844 interactions observed by
mass-spectrometry based analysis of protein complexes. It is updated monthly
and has been subsumed as part of the BioGRID project which, in turn is a broad
interaction database including more than 1.7m individual interactions [10].

Both these contain kinase (and phosphatase) target associations with the
latter (Yeast Kinome db) also containing these at the phosphosite level. As
will be shown below however, there was an issue of overlap between our
experimental facts and background knowledge which limited the predictions that
could be made.

####  2.2.1 Yeast perturbation observation facts and background knowledge
sources

From the perturbation study chosen, a total of 685,358 facts were extracted
containing the effect of the aforementioned 110 kinase and phosphatase
deletions over 13,258 individual phosphosites. They were extracted from a .csv
file containing the log2fold change recorded for each phosphosite between wild
type and $\delta$gene strains of yeast (i.e. strains where selected genes have
been deleted). For most phosphosites there were two replicates reported. Not
all phosphosites had both replicates for all $\delta$gene combinations. Only
phosphosites containing two replicates were considered and the mean value
between the two was taken for those that did. The code determined whether the
average of the log2 fold change values is less than or equal to -0.5, greater
than or equal to 0.5, or between those two values. Depending on the value of
the average, the code sets the value of a variable called “ change” to “down”,
“up”, or “unaffected”, respectively. Facts took the following form:

perturbs(’YCK2’, ’ENO2_pT324’, ’ENO2’, unaffected, 0.0475).

The background knowledge base includes facts that describe the relationships
between kinases, proteins and phosphosites at both the kinase-protein and
kinase-phosphosite level. These were taken from the two databases described
above and were transcribed into facts in the following format:

Protein level information: knownKtarget(’ARK1’, ’PAN1’).

Phosphosite level information: knowntarget(’ATG1’,’ALY1_pS813’).

The first set of facts describing kinase - protein level information came from
the Kinome Interaction Database and the Yeast KINOME database, while the
kinase - phosphosite level information came from the KINOME/bioGRID database
alone.

All scripts used during data and background knowledge curation and
transcription as well as .pl Prolog files can be found at
https://github.com/Dudelder/Yeast_Prolog_Compensatory). Total number of facts
that make up knowledge base can be seen in Table 1.

Facts List | Source |  | Size
---
(Number of facts)
| Dgene-Dgene phosphosite
---
phosphosite observations
| Perturbation
---
Study
685,358
Kinase Substrate Associations |  | Yeast Kinome
---
Database / bioGRID
7589
Kinase Protein Associations |  | Yeast Kinome
---
Database / bioGRID
1558
Kinase Protein Associations |  | Kinase Interaction
---
Database (lenient)
853
Kinase Protein Associations |  | Kinase Interaction
---
Database (Strict)
517
| Total = | 695,875
Table 1: Table listing the title, sources and associated number of facts that
make up the Background Knowledge section of the yeast logic model. The total
reported is for the strict cut-off of the Kinome Interaction Database. For the
lenient cut-off, the total is 696,211. Specific sources referenced [12], [8],
[2].

####  2.2.2 Overlap between experimental observations and background knowledge

There is an overlap between the observational facts at the kinase-phosphosite
level and the phosphosites in the background knowledge base. Specifically, out
of the total number of individual phosphosites, only 667 are found in the
kinase-phosphosite associations knowledge base (Figure 2 A), accounting for
approximately 5$\%$ of the total.. At the protein level things improve
slightly depending on the data source and cutoff used (for Kinase Interaction
Database). In the Kinome database (see Figure 2), there is an overlap between
the information obtained from Kinase Protein pairs and Kinase Protein
Phosphosite level triples. However, there is a discrepancy in the numbers. The
overlap consists of 72 deleted kinases when considering the Kinase Protein
pairs, but only 58 deleted kinases when extracting the information from the
Kinase Protein Phosphosite level triples. This difference may be due to errors
in formatting or the omission of a subset of data in one of the database
files. Looking at the Kinase Interaction database as a source, the overlap
depends on the cutoff used. Cutoffs of 6.73 and 4.72 correspond to a P value
of less than 0.01 and 0.05, respectively, for strict and lenient lists of
kinase-protein pairs as were suggested in the article. Accordingly the overlap
between the deleted kinases and the two lists of Kinase and their target
Protein pairs was 52 and 58 for the strict and lenient cutoffs, respectively.
As will be demonstrated in the results section of this paper, the poor overlap
was addressed by considering combinations of background knowledge sources.

Figure 2: Venn Diagrams depicting the overlap between our experimental data
and background knowledge sourced from Kinome Database / bioGRID [2]. In A)
total number of phosphosites is 13,256. B: Overlap between kinases in
KinomeDB/bioGRID and deleted kinases at the kinase-protein level of
association. C: Overlap between kinases in KinomeDB/bioGRID and deleted
kinases at the kinase-phosphosite level of association.

##  3 Results

###  3.1 Silence-of-K Rule iterations and Kinase Compensation Rule

Collecting the results was accomplished following a pipeline briefly described
here. The Prolog interpreter was used to perform queries across all kinase
deletions, and the resulting output was stored in lists. These lists were then
parsed by Python scripts to perform additional visualisation and analysis in
RStudio. This centered around the number of deletions that were predicted to
have worked. Additionally, RStudio was used to create the kinase compensation
networks, as described below in this section.

Using kinase-phosphosite level associations, the first iteration of the rule
proved to exhibit the best recall, confirming 21 kinase gene deletions having
taken place, thus confirming the ground truth established by PCR, as reported
in the study. Successive iterations of the rule, 2nd (checking uniqueness) and
3rd (checking majority of known phosphosites) were less accurate, both
confirming only 5 kinases as being silenced. Using kinase-protein level
associations the results improved overall depending on the source of kinase
protein pairs used. For the Kinome Database from bioGrid, 33, 7 and 3 kinase
deletions were confirmed for rule iterations 1, 2 and 3, respectively. Kinase
Interaction Database (KID) with a lenient cut-off yielded 26, 9 and 5, whereas
with a strict cut-off the results were 21, 8 and 4. Results summarised in
Table 2.

As the overlap with a ground truth (the deletion of a kinase encoding gene)
was poor, it was imperative to look at ways of improving this. Firstly,
combining the outputs of the two levels of background knowledge considered was
the most straightforward approach. Starting with the Kinome Database as a
Kinase-Protein association source, 54, 12 and 8 ground truths were captured
with each successive iteration of the rule. From the Kinase Interaction
Database, the lenient cut off gave 47, 14 and 10, and with the strict cut off
42, 13 and 9. Overall we noted a significant improvement which was mainly
evident in the first iteration of the rule, picking up increasingly more of
the ground truth. These can be considered True Positive results.

| Rule Iteration &
---
Background Knowledge
Level (Source where appropriate)
True Positive | False Negative
First Iteration & Phosphosite level | 21 | 63
Second Iteration & Phosphosite level | 5 | 79
Third Iteration & Phosphosite level | 5 | 79
First Iteration & Protein level (bioGrid) | 33 | 51
Second Iteration & Protein level (bioGrid) | 7 | 77
Third Iteration & Protein level (bioGrid) | 3 | 81
First Iteration & Protein level (KID Strict) | 21 | 63
Second Iteration & Protein level (KID Strict) | 8 | 76
Third Iteration & Protein level (KID Stict) | 4 | 80
First Iteration & Protein level (KID Lenient) | 26 | 68
Second Iteration & Protein level (KID Lenient) | 9 | 75
Third Iteration & Protein level (KID Lenient) | 5 | 89
Table 2: Table listing the True Positive and False Negative counts for each of
the rule iterations. The False Negatives refer to kinase gene deletions we
were not able to pick up due to either poor background information overlap or
rule iterations not performing as intended.

The next step was to combine all of the known background sources whilst
removing duplicates to avoid clashes. This yielded a kinase-protein pair
background knowledge base which covered (with at least one association) 79/82
of the deleted kinases in our dataset. This was the largest overall coverage
and yielded the best result, confirming 62 individual deletions as seen in
Figure 3. In contrast, iterations two and three yielded worse results, 13 and
8, respectively. From the above the best combinations between knowledge base
sources (KID: strict/lenient and KINOME:bioGRID), level of background
knowledge considered (kinase-phosphosite, kinase-protein or both) for each
rule iteration were:

  * •

Rule 1: All BK, both levels: 62 TP / 17 FN

  * •

Rule 2: KID Lenient cut-off, both levels: 14 TP / 46 FN

  * •

Rule 3: KID Lenient cut-off, both levels: 10 TP / 50 FN

Figure 3: Figure containing Venn Diagrams depicting the ground truth (true
positives) and false negatives captured by the top performing rule iterations
and background knowledge levels considered. In all three Venn diagrams (A,B
and C) the colour scheme is as follows: Purple for True Positives, light
yellow for False Negatives, light blue for background knowledge kinase protein
target pairs and light pink for gene deletion kinases. A: Rule iteration 1
with background knowledge at the protein level and from all sources combined
had the best recall with 62 True Positives identified out of a potential 78.
B: Rule iteration 2, Kinase Interaction database lenient cut-off and
phosphosite level information with 14 True Positives and 46 False Negatives
out of 60 potential. C: The third iteration of the rule performed best when
the background knowledge was sourced from KID with lenient cutoff at 10 True
Positives, 50 False Negatives out of a potential 64. The difference in
potential totals comes can be attributed to how many of the deleted kinases
were present and therefore queriable in each instance of background knowledge
source and protein or phosphosite level consideration.

####  3.1.1 Explaining the 17 missed deletions with a new rule

From the above it was decided to further focus on explaining the results
obtained from rule iteration 1, with the results of both levels included and
all background knowledge sources combined. However, as can be seen in the
triple Venn diagram of Figure 3 A, there were 17 kinase deletions the rule
iteration was not able to identify. As this iteration did not exclude shared
targets between kinases a possible hypotheses is that other kinases are acting
in a compensatory manner. For this the ’ykinasecompcheck/2’ rule was
developed.

In brief, it checks whether a kinase (K2) increases its activity when a
different kinase (DelK1) is silenced, in order to compensate for the loss of
K1 activity and maintain the overall integrity of the signalling pathway.
Specifically the Prolog notation the rule took the following form:

⬇

1ykinasecompcheck(DelK1, K2) :-

2 (perturbs(DelK1,_Pst,Tprot1,down,_) ;

3 perturbs(DelK1,_Pst,Tprot1,unaffected,_)),

4 sharedKinasetarget(DelK1,K2,Tprot1),

5 knowntarget(K2, Pst1, Tprot1),

6 (perturbs(DelK2,Pst1,Tprot1,up,_),

7 (DelK2 \= K2, DelK2 \= DelK1)).

Specifically, the first goal ’perturbs/5’ holds if the first argument (DelK1)
perturbs the fourth argument (TProt1) in the direction specified by the fifth
argument (down or unaffected). ’sharedKinasetarget/3’ is a goal that is
satisfied if the first argument (DelK1) and the second argument (K2) share the
third argument (Tprot1) as a target. The information used for this, in order
to maintain consistent background knowledge was the combination of all kinase
protein background knowledge sources. ’knowntarget/3’ is a goal that is
satisfied if Tprot1 is a known target of the K2 at the phosphosite specified
by the second argument (Pst1). Finally, (6) holds if the first argument
(DelK2) perturbs the second argument (Pst1) in the direction specified by the
fifth argument (up). The condition is true if DelK2 is different from both K2
and DelK1. The last goal ensures that kinase K2 indeed phosphorylates
phosphosite Pst1 (in the condition it, itself is not deleted).

Using the above rule and the list of 17 missed deleted Kinases, for each, a
further list of compensatory kinases was collated. These lists where then used
to build compensatory kinase interaction networks in order to visualise the
effect of kinase deletions on known kinase interaction networks. Predicted or
compensatory pathways were compared with pathways reconstructed from
literature [12]. Visualisations were carried out in Rstudio using the igraph
package.

None of the false negatives (kinase deletions not present in the rule
predictions) were present in the pathways associated with endocytosis,
transcription response and mating. Cell cycle had the largest number of
kinases compensating for ones from within our list. In Figure 4 Cell Cycle
pathway is highlighted with a red border, the High-Osmorality glucose (HOG)
pathway with a blue border and the meiosis with a fuchsia border. They are all
connected canonically, showing the relative accuracy of the rule described
above in picking up the shift in pathway activity to account for the deletion
of a given kinase within the pathway. For the highlighted pathways the
majority of their constituent kinases are present, deletions of which the
first iteration of the rule did not pick up. This indirectly suggests that in
response to certain key deletions these highlighted pathways are able to
rewire themselves in an effort to counteract the perturbation.

Figure 4: Predicted manner of kinase network rewiring following gene
deletions. Kinases within pathways act in a compensatory manner. A) Cell Cycle
pathway, predicted pathway rewiring. B) High-osmorality glucose and cell cycle
pathway rewiring. C) High-osmorality and Meiosis pathway rewiring. Highlighted
are the kinases belonging to the respective pathways that include False
Negative i.e. missed gene deletions within our predictions.

##  4 Discussion

Despite Saccharomyces cerevisiae being extensively studied as a model system,
there was a notable lack of coverage between kinase-phosphosite associations
and the observational (phosphoproteomics) data obtained from the gene knockout
study. Therefore implementing the methodology developed initially on the Toy
Model, proved to be more challenging than anticipated.

The approach that yielded the most accurate results, in terms of recall, with
the simplest iteration of the rule. The issue of poor overlap was also
partially overcome by combining background knowledge sources and levels of
granularity. Increasing complexity of rules did not yield better results.

Delving deeper into the False Negatives, i.e. kinases from verified deleted
genes which we were unable to pick up, yielded interesting results. The
interconnected nature of kinase interaction pathways was captured by the
’ykinasecompcheck/2’ rule. Using the 17 kinases as a base we were able to
identify and contextualise the majority of major large Yeast pathways, namely
HOG and cell cycle as well as smaller one such as Meiosis. For the remaining
pathways, their constituents were correctly captured as having been silenced.
Combining the outputs of the individual silenced gene queries offers an
insight into potential yeast signalling rewiring following strong perturbation
such as silencing of a group of kinases.

The issue with overlap can be addressed by initiatives such as bioGRID and
other datasets providing consistent information at both the protein and
phosphosite level interactions. This problem is even greater and more present
in the human kinome and proteome. However, combining resources, a crude
approach at best, seemed to alleviate this, at least partially. A recently
published study [6] made a significant contribution in assigning substrates to
a large number of active Ser/Thr kinases (approximately 84$\%$) in humans.

Next steps in applying this methodology would involve augmenting facts with
probabilities, via the ProbLog2 [4] framework. Specifically, the immediate aim
of this would be to capture the epistemic uncertainty in the kinase
phosphosite associations as well as the mix of alleatoric and epistemic
uncertainty that characterises the process of mass spectrometry based data
acquisition and related post processing. The resulting output could further
increase confidence in the predictions. Both of these require the mining of
appropriate and well curated information in order to prove viable.

## References

  * [1]
  * [2] Ashton Breitkreutz, Hyungwon Choi, Jeffrey R. Sharom, Lorrie Boucher, Victor Neduva, Brett Larsen, Zhen Yuan Lin, Bobby Joe Breitkreutz, Chris Stark, Guomin Liu, Jessica Ahn, Danielle Dewar-Darch, Teresa Reguly, Xiaojing Tang, Ricardo Almeida, Zhaohui Steve Qin, Tony Pawson, Anne Claude Gingras, Alexey I. Nesvizhskii & Mike Tyers (2010): _A global protein kinase and phosphatase interaction network in yeast_.  Science (New York, N.Y.) 328(5981), 10.1126/SCIENCE.1176495.  Available at https://pubmed.ncbi.nlm.nih.gov/20489023/.
  * [3] Amanda Clare & Ross D. King (2003): _Data Mining the Yeast Genome in a Lazy Functional Language_ , pp. 19–36.  10.1007/3-540-36388-2_4.  Available at http://www.aber.ac.uk/compsci/Research/bio/dss/polyfarm/.
  * [4] Anton Dries, Angelika Kimmig, Wannes Meert, Joris Renkens, Guy Van Den Broeck, Jonas Vlasselaer & Luc De Raedt (2015): _ProbLog2: Probabilistic logic programming_.  10.1007/978-3-319-23461-8_37.  Available at https://dtai.cs.kuleuven.be/problog.
  * [5] Antonio Goncalves, Irene Ong, Jeffrey A. Lewis & Vitor Santos Costa (2014): _Towards using probabilities and logic to model regulatory networks_.  Proceedings - IEEE Symposium on Computer-Based Medical Systems, pp. 239–242, 10.1109/CBMS.2014.9.  Available at https://ieeexplore.ieee.org/document/6881883.
  * [6] Jared L. Johnson, Tomer M. Yaron, Emily M. Huntsman, Alexander Kerelsky, Junho Song, Amit Regev, Ting Yu Lin, Katarina Liberatore, Daniel M. Cizin, Benjamin M. Cohen, Neil Vasan, Yilun Ma, Konstantin Krismer, Jaylissa Torres Robles, Bert van de Kooij, Anne E. van Vlimmeren, Nicole Andrée-Busch, Norbert F. Käufer, Maxim V. Dorovkov, Alexey G. Ryazanov, Yuichiro Takagi, Edward R. Kastenhuber, Marcus D. Goncalves, Benjamin D. Hopkins, Olivier Elemento, Dylan J. Taatjes, Alexandre Maucuer, Akio Yamashita, Alexei Degterev, Mohamed Uduman, Jingyi Lu, Sean D. Landry, Bin Zhang, Ian Cossentino, Rune Linding, John Blenis, Peter V. Hornbeck, Benjamin E. Turk, Michael B. Yaffe & Lewis C. Cantley (2023): _An atlas of substrate specificities for the human serine/threonine kinome_.  Nature 613(7945):4, 10.1038/s41586-022-05575-3.  arXiv:https://arxiv.org/abs/2042.12345.  Note.
  * [7] Ross D. King, Kenneth E. Whelan, Ffion M. Jones, Philip G.K. Reiser, Christopher H. Bryant, Stephen H. Muggleton, Douglas B. Kell & Stephen G. Oliver (2004): _Functional genomic hypothesis generation and experimentation by a robot scientist_.  Nature 427(6971), pp. 247–252, 10.1038/nature02236.  Available at https://www.nature.com/articles/nature02236.
  * [8] Jiaming Li, Joao A. Paulo, David P. Nusinow, Edward L. Huttlin & Steven P. Gygi (2019): _Investigation of Proteomic and Phosphoproteomic Responses to Signaling Network Perturbations Reveals Functional Pathway Organizations in Yeast_.  Cell reports 29(7), p. 2092, 10.1016/J.CELREP.2019.10.034.
  * [9] G. Manning, D. B. Whyte, R. Martinez, T. Hunter & S. Sudarsanam (2002): _The protein kinase complement of the human genome_.  Science (New York, N.Y.) 298(5600), pp. 1912–1934, 10.1126/SCIENCE.1075762.  Available at https://pubmed.ncbi.nlm.nih.gov/12471243/.
  * [10] Rose Oughtred, Jennifer Rust, Christie Chang, Bobby Joe Breitkreutz, Chris Stark, Andrew Willems, Lorrie Boucher, Genie Leung, Nadine Kolas, Frederick Zhang, Sonam Dolma, Jasmin Coulombe-Huntington, Andrew Chatr-aryamontri, Kara Dolinski & Mike Tyers (2021): _The BioGRID database: A comprehensive biomedical resource of curated protein, genetic, and chemical interactions_.  Protein science : a publication of the Protein Society 30(1), 10.1002/PRO.3978.  Available at https://pubmed.ncbi.nlm.nih.gov/33070389/.
  * [11] Jason Ptacek, Geeta Devgan, Gregory Michaud, Heng Zhu, Xiaowei Zhu, Joseph Fasolo, Hong Guo, Ghil Jona, Ashton Breitkreutz, Richelle Sopko, Rhonda R. McCartney, Martin C. Schmidt, Najma Rachidi, Soo Jung Lee, Angie S. Mah, Lihao Meng, Michael J.R. Stark, David F. Stern, Claudio De Virgilio, Mike Tyers, Brenda Andrews, Mark Gerstein, Barry Schweitzer, Paul F. Predki & Michael Snyder (2005): _Global analysis of protein phosphorylation in yeast_.  Nature 2005 438:7068 438(7068), 10.1038/nature04187.
  * [12] Sara Sharifpoor, Alex N. Nguyen Ba, Ji Young Young, Dewald van Dyk, Helena Friesen, Alison C. Douglas, Christoph F. Kurat, Yolanda T. Chong, Karen Founk, Alan M. Moses & Brenda J. Andrews (2011): _A quantitative literature-curated gold standard for kinase-substrate pairs_.  Genome Biology 12(4), 10.1186/GB-2011-12-4-R39.  Available at https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3218865/.
  * [13] K. E. Whelan & R. D. King (2008): _Using a logical model to predict the growth of yeast_.  BMC Bioinformatics 9(1), pp. 1–16, 10.1186/1471-2105-9-97/TABLES/12.  Available at https://link.springer.com/articles/10.1186/1471-2105-9-97https://link.springer.com/article/10.1186/1471-2105-9-97.