Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 2005-0461(Print)
ISSN : 2287-7975(Online)
Journal of Society of Korea Industrial and Systems Engineering Vol.44 No.4 pp.154-168
DOI : https://doi.org/10.11627/jksie.2021.44.4.154

The Analysis of COVID-19 Pooled-Testing Systems with False Negatives Using a Queueing Model

Kilhwan Kim†
Department of Management Engineering, Sangmyung University
Corresponding Author : khkim@smu.ac.kr
17/11/2021 02/12/2021 06/12/2021

Abstract


COVID-19 has been spreading all around the world, and threatening global health. In this situation, identifying and isolating infected individuals rapidly has been one of the most important measures to contain the epidemic. However, the standard diagnosis procedure with RT-PCR (Reverse Transcriptase Polymerase Chain Reaction) is costly and time-consuming. For this reason, pooled testing for COVID-19 has been proposed from the early stage of the COVID-19 pandemic to reduce the cost and time of identifying the COVID-19 infection. For pooled testing, how many samples are tested in group is the most significant factor to the performance of the test system. When the arrivals of test requirements and the test time are stochastic, batch-service queueing models have been utilized for the analysis of pooled-testing systems. However, most of them do not consider the false-negative test results of pooled testing in their performance analysis. For the COVID-19 RT-PCR test, there is a small but certain possibility of false-negative test results, and the group-test size affects not only the time and cost of pooled testing, but also the false-negative rate of pooled testing, which is a significant concern to public health authorities. In this study, we analyze the performance of COVID-19 pooled-testing systems with false-negative test results. To do this, we first formulate the COVID-19 pooled-testing systems with false negatives as a batch-service queuing model, and then obtain the performance measures such as the expected number of test requirements in the system, the expected number of RP-PCR tests for a test sample, the false-negative group-test rate, and the total cost per unit time, using the queueing analysis. We also present a numerical example to demonstrate the applicability of our analysis, and draw a couple of implications for COVID-19 pooled testing.



대기행렬을 이용한 위음성률이 있는 코로나 취합검사 시스템의 분석

김 길환†
상명대학교 경영공학과

초록


    1. 서 론

    코로나19 바이러스에 의한 감염증이 전 세계적으로 유행 하는 현재, 코로나19 감염자를 빠르게 찾아내어 격리, 치료하 는 것은 방역의 매우 중요한 요소 중 하나이다. 코로나19 감염자를 빠르게 찾아내기 위해서는 코로나19 감염자를 신속하고 정확하게 진단하는 검사 시스템의 구축이 필요하 다. 코로나19의 진단은 항원 및 항체 반응 검사, 핵산 증폭 검사, 배양 검사 등이 있다[26]. 항원 및 항체 반응 검사는 5~15분 이내에 결과를 얻을 수 있는 검사법이지만 실제 감염자를 양성으로 진단하는 민감도(sensitivity)가 낮아 제한 적인 상황에서만 사용되고 있다. 배양 검사는 바이러스 배양 과 검사에 5~9일 이상의 시간이 소요되므로 감염자를 조기에 진단하기 어렵기 때문에 이 방법도 현재 코로나19 진단에 사용되지 않고 있다. 현재 코로나19 검사에 주로 사용되는 방법은 핵산 증폭 검사로서, 세계보건기구(WHO)에서 권장 하고 있는 실시간 역전사 중합효소 연쇄반응(rt RT-PCR; real-time reverse transcriptase polymerase chain reaction) 검사 법을 사용하여 진단을 하고 있다. 이 방법의 민감도는 98.2% 이상이며, 비감염자를 음성을 진단하는 특이도(specificity) 도 100%로 매우 정확한 검사 방법이다. 그러나 이 검사법은 3~6 시간이 소요되어 항원 및 항체 반응 검사법에 비해 진단 시간이 긴 단점이 있다. 또한 전 세계적으로 보면 RT-PCR 검사에 필요한 검사 키트와 시약의 공급 부족이 큰 이슈로 대두되고 있다. 이런 이유로, 검사 장비와 인력의 제약 하에서 폭증하는 코로나19 검사량에 대응하기 위하여, 검사의 정확도를 약간 희생하면서 여러 검체를 한 번에 검사하는 취합검사법(pooled test)을 RT-PCR에 적용하는 연구가 수행되었다[1, 9, 13, 16, 28].

    취합검사법은 도프만이 2차 세계 대전 중에 군 징집자 중 매독 감염자를 효율적으로 검출해내기 위한 방법으로 처음 제안한 검사 방법이다[10]. 이 방법에서는 일정 수 의 혈액 표본을 취합하여 그룹 검사한 후, 음성인 경우 취합된 혈액을 모두 음성으로, 양성인 경우 취합된 혈액 표본을 개별 검사하여 양성자를 찾아낸다. 검사 대상의 확진율(검사 대상 대비 확진자의 비율)이 매우 낮은 경 우에는 취합된 표본이 크더라도 그룹 검사에서 음성일 가능성이 커져서 취합검사법을 사용하면 더 적은 검사 횟수로 검사를 완료할 수 있다. 그러나 확진율이 높은 경 우에는 취합된 표본이 크면 그룹 검사에서 양성일 확률 이 커져서 취합검사 후에 수행해야 하는 개별 검사로 인 해 오히려 검사 횟수가 증가하게 된다. 그렇기 때문에 검 사 대상의 확진율에 따라 그룹 검사를 하는 검체의 수를 최적으로 설정하는 것이 매우 중요하다. 도프만의 연구 이후, 다양한 취합검사법이 제안되었다. 대표적으로 그 룹 검사에서 양성이 나온 경우 개별 검사를 바로 하는 것이 아니라 소규모 그룹으로 다시 나누어 검사한 후, 소 규모 그룹에서 양성이 나온 그룹만 다시 개별 검사를 하 는 방법 등도 제안되었다. 다양한 취합검사법에 대해서 는 Du[12]를 참고하기 바란다. 도프만의 연구 이 후, 취 합검사는 매독 이외에도 HIV나 다른 바이러스와 박테리 아의 검출에도 널리 사용되고 있다[19, 24, 27].

    하지만 취합검사법에 대한 대부분의 연구는 검사 대 상이 모두 확정되어 있는 상황에서 검사의 정확도와 검 사 시스템의 효율적 사용의 관점에서 최적의 취합 표본 수를 결정하는 것에 주안점을 두었다. 도프만이 제안했 던 군 징집자의 매독 검사처럼 검사 대상의 수가 어느 정도 확정되어 있는 경우에는 이러한 방법이 적용될 수 있으나, 검사 대상의 발생이 확률적으로 발생하는 경우 에는 검사 요구 발생의 확률적 속성과 검사 시간의 확률 적 변동성을 고려하여 최적 취합 표본 수를 결정하는 것 이 필요하다. 특히, 코로나19처럼 정확한 진단만큼이나 신속한 진단이 필요한 검사에서는 검사의 정확도와 검사 시스템의 효율적 운용뿐 아니라, 검사 요청이 의뢰된 시 점부터 검사 결과가 통보될 때까지의 처리 시간도 검사 시스템의 중요한 성과 지표로 고려되어야 한다.

    검사 요청이 확률적으로 발생하고 검사 요청의 처리 시간의 단축이 중요한 경우, 확정적인 상황을 가정한 최 적 취합검사 수는 더 이상 최적이 아닐 수 있다.예를 들 어 매우 낮은 양성률을 보이는 경우, 확정적인 상황에서 는 취합검사 수를 크게 설정하는 게 최적이 된다. 그러나 검사 요청이 확률적으로 발생하고 검사 요청의 발생 주 기가 긴 경우에는, 이미 도착한 검사 요청들이 취합검사 에 필요한 검체 수에 도달될 때까지 검사를 기다려야 하 므로 처리시간의 평균과 분산이 커지게 된다. 코로나19 처럼 강한 전염성이 있는 경우, 확진의 지연은 전파의 가 능성을 크게 키우고 자가 격리 대상자의 불편을 증가시 키게 된다. 따라서 검사 요청의 발생과 검사 시간의 불확 실성 하에서 검사 요청의 처리 시간을 분석하여 최적 취 합 표본 수를 결정하는 방법이 필요하다.

    이러한 이유로, Abolnikov and Dukhovny[2]는 HIV 바 이러스 취합검사 시스템에 대하여 검사 요청의 도착을 포아송 과정으로 가정하여, 배치의 크기에 따라 서비스 시간이 달라지는 배치 서비스 대기행렬 모형을 사용하 여 검사 요청의 처리시간을 분석하였다. Abolnikov and Dukhovny[2]의 연구 이후, Bar-Lev et al.[4]와 Claeys et al.[8] 등도 배치의 크기에 따라 서비스 시간이 달라지는 배치 서비스 대기행렬을 사용하여 확률적 상황 하의 취 합검사 시스템의 성능 분석을 수행하였다. 특히, 이 연구 에서는 검사 요청의 도착 과정과 확진율을 고려하여 최 적의 최소 취합검사 수 a와 최대 취합검사 수 b를 결정 하는 방법을 제시하였는데, Bar-Lev et al.[4]는 검사 요 청이 개별적으로 도착하는 포아송 과정을 따른다고 가정 하여 M/G(a,b)/1 대기행렬로, Claeys et al.[8]은 취합검사 와 개별 검사 시간이 확정적이라고 가정하고 이산 시간 모형인 GeoX/G(a,b)/1 대기행렬 모형으로 최적 취합검사 수 ab를 탐색하였다. 아울러 직접적으로 취합검사 시 스템을 목표로한 연구는 아니지만, 취합검사 시스템의 분석에 활용할 수 있는, 배치의 크기에 따라 서비스 시간 이 달라지는 배치 서비스 대기행렬에 대한 연구도 널리 수행되어 왔다[3, 5-7, 11, 20, 22, 23, 25, 30].

    그러나 지금까지의 대기행렬을 이용한 취합검사 시스 템에 대한 분석은 모두 확률적 도착 과정과 확진율만을 고려하고, 검사의 민감도와 특이도는 취합검사 수에 무 관하게 항상 정확하다고 가정을 하였다[2, 4, 8]. 즉, 감 염자이지만 음성으로 판단되는 위음성률(false negative rate)이나, 비감염자이지만 양성으로 판단되는 위양성률 (false positive rate)이 0이라고 가정되었다. 그러나 코로 나19 rt RT-PCR 검사는 민감도가 항상 100%가 되기 어 려우며, 최근의 코로나19 바이러스에 대한 취합검사법에 대한 연구에서도 알 수 있듯이 취합되는 검체의 수가 증 가할수록 검출해야 하는 유전자가 희석되어 검사의 민감 도가 저하되는 현상이 발생한다[1, 9, 13, 16, 28]. 따라서 취합검사 시스템의 성능을 분석할 때, 검사 요청의 처리 시간뿐 아니라, 취합검사 수에 따른 검사의 정확도의 변 화를 모두 고려하는 것이 바람직하다.

    본 연구는 검사 요청의 도착이 확률적으로 발생하는 상황에서 취합검사 시스템의 성능을 분석하기 위하여 기 존의 연구처럼 배치의 크기에 따라 서비스 시간이 달라 지는 배치 서비스 대기행렬을 이용하지만, 취합검사 수 에 따라 검사의 민감도가 달라지는 점을 모형화하여 최 적 취합검사 수 (a,b)를 분석한다. 따라서 본 논문은 취 합검사에 따른 민감도의 저하-위음성률의 증가-를 고려 하지만 검사 요청이 이미 확정적으로 정해져 있는 상황 을 가정하여 최적 취합검사 수를 탐색한 연구나[1, 9, 13, 28], 검사 요청의 확률적 도착 과정을 고려하였지만 취합 검사에 따른 민감도의 저하를 고려하지 않고 항상 검사 의 민감도가 100%라고 가정한 연구[2, 4, 8], 그리고 시 뮬레이션 기법으로 취합검사 시스템을 탐색한 연구 [16] 등과 차별성을 가진다고 할 수 있다.

    이 논문의 나머지 부분은 다음과 같이 구성되어 있다. 제2장에서는 위음성률(false negative rate)이 있는-민감도 가 100%가 아닌-코로나 취합검사 시스템을 대기행렬로 모형화한다. 제3장에서는 대기행렬 이론을 이용하여 임 의 시점에 검사 시스템에 있는 검사 요청 수의 분포를 도출한다. 제4장에서는 취합검사 시스템의 성능 지표를 도출한다. 제5장에서는 수치 예제를 이용하여 다양한 확 진율과 검사 요구의 도착률에 대해 취합검사 수에 따라 취합검사 시스템의 성능이 어떻게 변하는지를 탐색해 본 다. 제6장에서는 연구의 결과의 의미를 토의하고 향후 연구 방향을 모색해 본다.

    2. 코로나 취합검사 시스템의 대기행렬 모형

    이 논문에서는 다음과 같은 코로나 취합검사 시스템을 가정한다.

    취합검사 시스템에는 J 대의 유전자 증폭 검사 장비가 있다. 검사가 필요한 검체는 도착률 λJ의 포아송 과정으로 도착한다. 도착한 검체는 J 대의 검사 장비 중 무작위로 추첨된 장비로 전달된다. 따라서 검체는 λ = λJ/J의 도착 률을 가지는 독립적인 포아송 과정으로 각 검사 장비에 도착하게 된다. 따라서 앞으로 취합검사 시스템을 분석 할 때는 λ의 도착률로 검체가 도착하는 검사 장비 한대 의 성능을 분석할 것이다. 아울러 도착한 검체가 실제 양 성일 확률은 p로, 도착 과정이나 다른 검체가 양성일 확 률과 서로 독립이라고 가정한다.

    검사 장비에 도착한 검체는 검사 장비가 검사 중이면 대기열에서 대기한다. 대기열의 용량은 무한이라 가정한 다. 아울러 검사 장비가 유휴 중이라도 도착한 검체 수가 a개 미만이면 취합검사를 위해 대기한다. 단, a ≥ 1. 검 사장비가 유휴 중이거나 취합검사를 마친 시점에서, 시 스템에서 “대기 중인” 검체 수가 a개 이상이 되면 바로 취합검사를 수행하는데, 검사의 효율과 검사 장비의 용 량을 고려하여 취합검사 수는 최대 b개로 제한된다. 단, 1 ≤ ab. 따라서 취합검사에서 한 번에 취합하여 검사 하는 검체의 수를 Y 라고 하면, 검체 도착의 확률적 성질 때문에 Y는 확률 변수가 되고, aYb이다.

    Y개의 검체에 대한 취합검사는 먼저 Y 개의 검체를 하나의 그룹으로 묶어 그룹 검사를 실시한다. 그룹 검사 에서 음성으로 판정되면 취합검사에 포함된 모든 검체를 음성으로 판정한다. 코로나19 유전자 증폭 검사의 특이 도(음성을 음성으로 판정하는 비율)은 100%이기 때문에 그룹 검사의 위양성 진단은 없다고 가정한다.

    그룹 검사에서 양성으로 판정되면 Y개의 검체를 개별 검사한다. 문제의 단순화를 위하여 그룹 검사가 양성인 데 모든 개별 검사가 음성인 경우에는 개별 검사의 정확 성을 신뢰하여 재검사를 수행하지 않는다고 가정한다.

    2.1 그룹 검사의 민감도

    코로나19 유전자 증폭 검사의 민감도에 대한 확률적 속성을 다음처럼 모형화한다.

    RT-PCR 검사는 목표유전자를 검출하기 위하여 유전자 증폭을 반복하여 검출하고자 하는 유전자가 있는지를 검사 한다. 일반적으로 한 번의 유전자 증폭 사이클 동안 검출해야 할 목표유전자가 두 배로 증가한다. 목표유전자의 검출은 목표유전자와 결합하는 형광물질의 형광신호로 검출하는 데, 형광신호가 역치(threshold) 이상이 되는 사이클의 수 Ct(cycle threshold)를 측정하여 양성 여부를 결정한다 .일반 적으로 검사 방법에 따라 최대 증폭 사이클 수 cmax가 정해져 있어서, Ct 값이 cmax 이하이면 양성으로 확진하고, cmax 사이클만큼 유전자 증폭을 하였는데도 형광신호가 역치에 도달하지 않으면 음성으로 판정한다. 일반적으로 cmax는 35~40 정도로 설정된다[15, 21, 26, 29].

    Mutesa et al.[21]의 연구에 따르면, 그룹 검사에 포함 되는 양성 검체 수가 오직 하나이면, 그룹 검사에 포함 된 전체 검체 수 Y가 2 배 증가할 때마다 검출해야 할 목표유전자가 1/2로 희석되기 때문에 유전자 검출을 위 한 사이클 역치 Ct가 1씩 증가할 것으로 예측하였고, 실제 실험 결과에서도 목표유전자가 유전자 N인 경우 ΔCt = 1.0 ± 0.15(평균 1, 표준편자 0.15)만큼 증가하였고, 목표유전자가 orflab인 경우에는 ΔCt = 1.1 ± 0.14만큼 증 가하였음을 밝혔다. 마찬가지로 그룹 검사된 표본에 양 성 검체가 K개 있다고 하면, 목표 유전자의 농도가 K 배만큼 증가하므로 더 적은 유전자 증폭의 사이클 수에 서 형광 검출이 될 것이다. 따라서, Y 개의 검체를 취합 한 표본에 K개의 양성 검체가 있을 때의 사이클 역치 CTY/K는 다음과 같이 최종적으로 표현할 수 있다[16].

    C t Y / K = C t + log 2 ( Y/K ) · Δ Ct.

    이 논문에서는 Mutesa et al.[21]의 연구의 가정에 따라, 개별 양성 검체의 사이클 역치 Ct와, 그룹 검사의 검체 수가 2배로 늘 때마다 늘어나는 사이클 역치 증가량 ΔCt가 모두 독립적인 정규 분포를 따른다고 가정한다.

    C t N ( μ , σ ) , Δ C t N ( Δ μ , Δ σ ) .

    그러므로 CTY/K는 다음과 같은 정규 분포를 따른다.

    C t Y / K N ( μ + log 2 ( Y K ) · Δ μ , σ + log 2 ( Y K ) · Δ σ ) .
    (1)

    단, Y ≥ 1, 1 ≤ KY .

    그룹 검사의 검체 수가 Y = y일 때, 한 그룹 검사 표 본에 포함된 양성 검체의 수 K는 확진율이 p일 때 다음 과 같은 이항 분포를 따른다.

    Pr { K = k | Y = y } = ( y k ) p k ( 1 p ) y k , 0 k y .

    그룹 검사의 검체 수가 Y = y일 때, 그룹 검사에서 음 성으로 판독하는 최대 증폭 사이클 수를 c max y 라고 하고, 그룹 검사의 검체 수가 Y = y일 때의 그룹 검사의 민감 도를 1 - αy라고 하자. 그러면 그룹 검사의 민감도는 양 성 검체가 포함되었을 때(K ≥ 1) Cty/K가 최대 증폭 사 이클 수 c max y 보다 작거나 같아서 양성으로 판독하는 확 률이므로 다음처럼 표현할 수 있다.

    1 α y = k = 1 y Pr { K = k | Y = y } · Pr { C t y / k c y max } k = 1 y Pr { K = k | Y = y } = k = 1 y ( y k ) p k ( 1 p ) y k · Φ ( c y max μ log 2 ( y/k ) · Δ μ σ + log 2 ( y/k ) · Δ σ ) 1 ( 1 p ) y
    (2)

    단, Φ(∙)는 표준정규분포의 CDF(Cumulative Distribution Function)이다.

    여기서 주의할 점은, y 개의 검체를 취합검사할 때의 개별 앙성 검체의 민감도는 1 - αy와 같지 않다는 것이 다. 왜냐하면 개별 양성 검체가 양성으로 최종 판별되기 위해서는 그룹 검사에서 양성으로 판별된 후, 개별 검사 에서 다시 양성으로 판별되어야 하기 때문이다. 본 논문 에서는 그룹 검사의 민감도 1 - αy가 특정 값 1 - αG이상 이 되도록 c max y 를 설정한다고 가정한다. 예를 들어 그룹 검사가 개별 검사의 민감도에 비해 99% 이하로 민감도 가 저하되는 것을 방지한다면 c max y 는 1 - αy가 0.99 이상 이 되는 가장 작은 자연수가 된다. 따라서 c max y 는 다음처 럼 결정된다.

    min c max y { 1 , 2 , } c max y s . t . α y α G .
    (3)

    2.2 취합검사 시간의 분포

    그룹 검사의 민감도에 대한 제약 조건 (3)이 만족되면 취합검사 시간의 분포는 다음처럼 모형화할 수 있다.

    y개의 검체를 취합하여 검사하는 취합검사의 총 검사 시간을 Sy라고 하자. 총 검사 시간 Syy개의 검체를 그룹 검사하는 시간 S y G 와 그룹 검사에서 양성일 경우 y 개의 개별 검체에 수행되는 총 개별 검사 시간 S y I 로 구성된다. 아울러 IG를 그룹 검사가 양성으로 판정되면 1, 음성으로 판정되면 0의 값을 갖는 지시 변수로 정의하자.

    그룹 검사나 개별 검사 모두 유전자 증폭 검사를 준비하 는 시간 Tr과 유전자 증폭 사이클을 반복하는 유전자 증폭 검사 시간으로 구성된다고 가정한다. N을 그룹 검사에서 수행한 유전자 증폭검사 사이클의 횟수라고 하고, Tc,nn 번째 유전자 증폭 사이클을 수행하는 데 걸린 시간이라고 하면, 그룹 검사 시간 S y G 는 다음처럼 표현된다.

    S y G = T r + T c , 1 + T c , 2 + + T c , N .

    유전자 증폭 검사를 준비하는 시간 Tr과 유전자 증폭 검사의 한 사이클을 수행하는데 걸리는 시간은 Tc,n은 검사의 종류(그룹 검사 또는 개별 검사)나 유전자 증폭 횟수 N에 무관하게 모두 독립인 다음과 같은 지수분포 를 따른다고 가정한다.

    T r exp ( μ r ) , T c , n exp ( μ c ) .

    만약 Y = y 개의 취합 검체에서 양성인 검체의 수가 K = 0이라면, N = c y max 번 사이클 동안 유전자 증폭 검사 후에 음성으로 판정될 것이다. 따라서 Y = y이고 K = 0 일 때의 S y G 의 조건부 LST(Laplace-Stieltjes Transform)와 IG의 조건부 확률은 다음처럼 표현된다.

    E [ e s S y G | K = 0 , Y = y ] = μ r μ r + s ( μ c μ c + s ) c max y Pr { I G = 0 | K = 0 , Y = y } = 1 Pr { I G = 1 | K = 0 , Y = y } = 1.

    만약 Y = y개의 취합 검체에서 양성인 검체의 수가 K = k ≥ 1이라면, 사이클 역치 Cty/k의 값에 따라 그룹 검사 시간 S y G 가 달라질 것이다. 만약 C t y / k > c y max 사이클 동안 유전자 증폭 검사 후에 음성으로 잘못 판정될 것이고, n - 1 < Cty/kn이면 (n = 1,2,…, c max y ), n 사이클 동안 유전자 증폭 검사한 후에 양성으로 판정될 것이다. Cty/k ≤ 0이면 유전자 증폭을 하기 전에 목표 유 전자가 검출된다고 가정한다. 따라서 Y = y이고 K = k일 때의 S y G IG의 결합분포는 다음처럼 표현된다.

    Pr { I G = 0 | K = k , Y = y } E [ e s S y G | K = k , Y = y ] = ( 1 γ y / k ) μ r μ r + s ( μ c μ c + s ) c max y Pr { I G = 1 | K = k , Y = y } E [ e s S y G | K = k , Y = y ] = μ r μ r + s n = 0 c max y γ y / k , n ( μ c μ c + s ) n .

    단, 1 ≤ ky이고,

    γ y / k = Pr { I G = 1 | K = k Y = y } = 1 Pr { I G = 0 | K = k , Y = y } = Φ ( c max y μ log 2 ( y/k ) · Δ μ σ + log 2 ( y/k ) · Δ σ ) γ y / k , 0 = Pr { I G = 1 , N = 0 | K = k Y = y } = Φ ( μ log 2 ( y/k ) · Δ μ σ + log 2 ( y/k ) · Δ σ ) γ y / k , n = Pr { I G = 1 , N = n | K = k Y = y } = Φ ( n μ log 2 ( y/k ) · Δ μ σ + log 2 ( y/k ) · Δ σ ) Φ ( n 1 μ log 2 ( y/k ) · Δ μ σ + log 2 ( y/k ) · Δ σ ) , 1 n c max y .

    따라서 Y = y이고 K = k일 때의 S y G 의 조건부 LST와 IG의 조건부 확률은 다음처럼 표현된다.

    E [ e s S y G | K = k , Y = y ] = μ r μ r + s ( μ c μ c + s ) c y max { ( 1 γ y / k ) + n = 0 c max y γ y / k , n ( μ c + s μ c ) c y max n } Pr { I G = 0 | K = k , Y = y } = 1 Pr { I G = 1 | K = k , Y = y } = 1 γ r / k .

    단, 1 ≤ ky.

    그룹 검사에서 양성으로 판별되면 개별 검사에 들어 간다. 개별 검체의 검사 시간도 유전자 증폭 검사를 준비 하는 시간 Tr과 유전자 증폭 사이클 시간으로 구성된다. 개별 검사의 준비 시간과 유전자 증폭 사이클 시간은 개 별 검사 간에 서로 독립이고, 그룹 검사 시간의 준비시간 과 유전자 증폭 시간과도 서로 독립이라고 가정한다. 한 검체의 개별 검사 시간은 하나의 검체로만 이루어진 그 룹검사 시간과 동일하므로 S 1 G 이 된다. 따라서 y 개의 취 합 검체 중에서 k 개의 양성 검체가 있다고 하면, 전체 개별 검사 시간 S y I k 개의 양성 검체에 대한 개별 검 사 시간과 y - k 개의 음성 검체에 대한 개별 검사 시간 의 합이 되므로, 양성 검체가 k 개 일 때의 총 개별 검사 시간 S y I 의 조건부 LST S y I * ( s ; k ) 는 다음처럼 표현된다.

    S y I * ( s ; k ) = E [ e s S y I | I G = 1 , K = k , Y = y ] = E [ e s S 1 G | K = 1 , Y = 1 ] k · E [ e s S 1 G | K = 0 , Y = 1 ] y k = ( μ r μ r + s ) y ( μ c μ c + s ) y c max × { ( 1 γ 1 / 1 ) + n = 0 c max γ 1 / 1 , n ( μ c + s μ c ) c max n } k .

    단, 1 ≤ ky. 위의 결과와 그룹 검사 시간 S y G 의 조건부 LST에서 Sy의 LST S y * (s) 를 다음처럼 구할 수 있다.

    S y * ( s ) = k = 0 y Pr { K = k | Y = y } × [ Pr { I y = 0 | K = k , Y = y } E [ e s S G | K = k , Y = y ] + Pr { I y = 1 | K = k , Y = y } E [ e s ( S G + S I ) | K = k , Y = y ] ] = μ r μ r + s ( μ c μ c + s ) c y max × [ ( 1 p ) y + k = 1 y ( y k ) p k ( 1 p ) y k × { ( 1 γ y / k ) + S y I * ( s ; k ) · n = 0 c y max γ y / k , n ( μ c + s μ c ) c y max n } ] .
    (4)

    (4)를 s에 대하여 미분하여 s = 0을 대입하면 다음을 얻는다.

    E [ S y ] = 1 μ r + c y max μ c + k = 1 y ( y k ) p k ( 1 p ) y k { γ y / k E [ S y I ] k n = 0 c max y γ y / k , n ( c y max n ) μ c }
    (5)

    E [ S y 2 ] = 2 μ r 2 + 2 c y max μ r μ c + c y max ( c y max + 1 ) μ c 2 + 2 ( 1 μ r + c y max μ c ) k = 1 y ( y k ) p k ( 1 p ) y k × { γ y / k · E [ S y I ] k n = 0 c max y γ y / k , n ( c y max n ) μ c } + k = 1 y ( y k ) p k ( 1 p ) y k { γ y / k · E [ S y I 2 ] k 2 E [ S y I ] k n = 0 c max y γ y / k , n ( c y max n ) μ c + n = 0 c max y γ y / k , n ( c y max n ) ( c y max n 1 ) μ c 2 } .
    (6)

    단,

    E [ S y I ] k = S y I * ( 0 , k ) = y ( 1 μ r + c max μ c ) k μ c n = 0 c max γ 1 / 1 , n ( c max n ) E [ S y I 2 ] k = S y I * ( 0 , k ) = y ( y + 1 ) μ r 2 + 2 y 2 c max μ r μ c + y c max ( y c max + 1 ) μ c 2 2 k y μ c ( 1 μ r + c max μ c ) n = 0 c max γ 1 / 1 , n ( c max n ) + k ( k 1 ) μ c 2 ( n = 0 c max γ 1 / 1 , n ( c max n ) ) 2 + k μ c 2 n = 0 c max γ 1 / 1 , n ( c max n ) ( c max n 1 ) .

    2.3 취합검사 당 검사 횟수의 분포

    Myy개의 검체로 취합 검사를 할 때 취합검사 1회 당 발생하는 검사 횟수라고 하자. 그룹 검사에서 음성으로 판정되면 개별 검사는 진행되지 않으므로 My = 1이 될 것이다. 반면 그룹 검사에서 양성으로 판정되면 취합된 검체 수만큼 의 개별 검사가 진행되므로 My = 1 + y가 될 것이다. 따라서 취합검사의 크기가 Y = y일 때, 취합된 검체에서 양성 검체 의 수가 K = k라면, 취합검사 1회 당 검사 횟수 My의 분포는 다음과 같다.

    Pr { M y = 1 | K = k } = 1 Pr { M y = 1 + y | K = k } = { 1 , k = 0 1 γ y / k , 1 k y .

    그러므로 취합검사의 크기가 Y = y일 때의 취합 검사 당 검사 횟수 My의 분포는 다음과 같다.

    Pr { M y = 1 } = 1 Pr { M y = 1 + y } = ( 1 p ) y + k = 1 y ( y k ) p k ( 1 p ) y k ( 1 γ y / k ) .

    따라서 My의 기대치는 다음과 같다.

    E [ M y ] = 1 + y k = 1 y ( y k ) p k ( 1 p ) y k γ y / k .
    (7)

    2.4 그룹 검사에서 위음성으로 진단되는 양성 검체 수의 분포

    코로나19 바이러스에 대한 유전자 증폭 검사의 특이 도는 100%이므로 위양성(음성을 양성으로 오진단)에 따 른 사회적 불편 비용은 없으나, 민감도는 98% 정도로 위 음성(양성을 음성으로 오진단)에 따른 사회적 위험 비용 이 발생할 수 있다. 취합검사에서 양성 검체가 위음성으 로 잘못 판단되는 사건은 그룹 검사에서 양성 검체가 있 음에도 음성으로 판정되는 경우와, 그룹 검사에서 양성 으로 판정되어 개별 검사를 하였지만 개별 검사에서 양 성 검체가 음성으로 오진단되는 경우가 있다. 개별 검사 의 위음성률은 취합검사를 수행하거나 하지 않거나 동일 한 비율로 발생하므로, 이 절에서는 개별 검사에서 위음 성으로 오진단되는 경우는 제외한, 취합검사를 도입함으 로써 그룹 검사에서 위음성으로 오진단되는 양성 검체 수의 분포에 대해서만 분석한다.

    y개의 취합 검체에 k개의 양성 검체가 있을 때 그 룹 검사에서 위음성으로 잘못 판정될 확률은 다음과 같다.

    Pr { I G = 0 | Y = y , K = k } = 1 γ y / k , 1 k y .

    따라서 y개의 검체를 취합하여 검사했을 때 그룹 검사에 서 위음성으로 진단되어 진단 체계를 벗어나는 양성 검체의 수를 Fy라고 하면, Fy의 분포는 다음과 같이 표현된다.

    Pr { F y = k } = Pr { K = k | Y = y } · Pr { I G = 0 | Y = y , K = k } = ( y k ) p k ( 1 p ) y k · ( 1 γ y / k ) .

    그러므로 y개의 검체를 취합하여 검사했을 때 그룹 검사에 의해 위음성으로 오진단되는 양성 검체 수 Fy의 기대치는 다음과 같다.

    E [ F y ] = k = 1 y k · Pr { K = k , I G = 0 | Y = y } = y p k = 1 y k · ( y k ) p k ( 1 p ) y k · γ y / k .
    (8)

    3. 취합검사 시스템 내의 검체 수 분포

    2절에서 기술된 취합검사 시스템의 검체 수 분포는 M/G(a,b)/1 대기행렬로 분석될 수 있다. M/G(a,b)/1 대기 행렬은 배치 서비스를 수행하는 대기행렬로, 배치 서비 스의 크기가 서비스 시작 시점의 고객수에 따라 a에서 b 까지 변화하고, 배치의 크기에 따라 서비스 시간이 변하 는 대기행렬이다. 단, 1 ≤ ab. 배치의 크기가 변화하 는 배치 서비스 대기행렬에 대한 연구는 여러 연구자에 의해 수행되어 왔다[3, 5-7, 11, 20, 22, 23, 25, 30]. 본 연 구에서는 기존 연구의 분석 방법을 이용하여 먼저 취합 검사 종료 직후 시점에 시스템에 있는 검체 수 분포를 구한 후, 이를 이용하여 임의시점에 시스템에 있는 검체 수 분포를 구하고자 한다.

    Lmm 번째 취합검사 종료 직후 시점(취합검사하 는 모든 검체에 대한 그룹검사와 필요한 개별검사를 모 두 완료한 시점)에 시스템에 남아 있는 검체 수라고 하 자. 단, m = 1,2,…. 그리고 Lm의 극한 확률(limiting probabilities)과 PGF(Probability Generating Function)를 π i = lim m P r { L m = i } Π ( z ) = i = 0 π i z i 로 정의하자.

    취합검사 종료 직후 시점에 시스템에 남아 있는 검체 수 Lm에 따라 다음 취합검사의 크기 Y가 결정된다. 만 약 취합검사 종료 직후 시점에 시스템에 남아 있는 검체 수가 a 미만이면, 검체 수가 a가 될 때까지 대기한 후, a 개의 검체가 취합검사된다. 이때의 취합검사 시간은 Sa가 된다. 만약 취합검사 종료 직후 시점에 검체가 y개이고 ayb - 1이면, y 개의 검체가 모두 취합검사되고 취 합검사 시간은 Sy가 된다. 만약 취합검사 종료 직후 시 점에 시스템에 검체 수가 b 이상이면 최대 그룹 검사 크 기인 b로 취합검사가 이루어지고, 그 때의 취합검사 시 간은 Sb이다. 그러므로 m-번째 취합검사 종료 직후 시점 과 다음 취합검사 종료 직후 시점 사이에 검사 받고 나 가는 검체 수와 새롭게 도착하는 검체 수를 고려하면, 다 음과 같은 관계식을 구할 수 있다.

    L m + 1 = { A a , 0 L m a 1 A L m , a L m b 1 L m b + A b , L m b .

    단, Ayy개의 검체를 취합하여 검사하는 시간 Sy 동 안 새롭게 시스템에 도착한 검체 수. 따라서 안정 상태에 서 취합검사 종료 직후 시점의 검체 수의 PGF Π(z)는 다음처럼 표현된다.

    Π ( z ) = i = 0 a 1 π i A a ( z ) + i = a b 1 π i A i ( z ) + i = b π i z i b A b ( z ) .

    단, A y ( z ) = E [ z A y ] = S y * ( λ λ z ) , ayb. 위의 식을 Π(z)에 대해여 정리하면 Π(z)에 대한 다음 식을 얻는다.

    Π ( z ) = i = 0 a 1 π i ( z i A b ( z ) z b A a ( z ) ) + i = a b 1 π i ( z i A b ( z ) z b A i ( z ) ) A b ( z ) z b .
    (9)

    (9)에서 λE [Sb] < b이면 검체 수에 대한 안정 상태 확 률이 존재하며, Π(z)의 분모식 Ab(z) = zb는 |z| ≤ 1 단위 원 안과 경계에 b 개의 해를 갖는다. 단위원 안과 경계에 있는 b 개의 해를 z0, z1, z2, ⋯, zb-1 이라고 하자. 그러면 Rouche의 정리에 의해 그중 하나는 z0 = 1이고 나머지 b - 1개의 해는 분자식에도 해가 되어야 한다[5]. 즉, j = 0,1,…,b - 1에 대해 Ab(zj) = z j b 이고 (9)의 분자식에서 도 zj가 해이므로 다음이 성립한다.

    i = 0 a 1 π i ( z j i A a ( z j ) ) + i = a b 1 π i ( z j i A i ( z j ) ) = 0 , j = 1 , 2 , , b 1.
    (10)

    그리고 확률의 정규 조건에서 lim z 1 Π ( z ) = 1 이어야 하 므로 다음 식을 얻는다.

    i = 0 b 1 π i i + ( 1 i = 0 b 1 π i ) b = i = 0 a 1 π i λ E [ S a ] + i = a b 1 π i λ E [ S i ] + ( 1 i = 0 b 1 π i ) λ E [ S b ] .
    (11)

    z1, z2, …, zb-1 이 서로 다른 해이면, (10)와 (11)에서 π0, π1, …, πb-1 에 대한 b개의 관계식을 얻을 수 있으므 로 (9)의 분자식의 미지수 π0, π1, …, πb-1 을 구할 수 있 다. 만약 z1, z2, …, zb-1 에서 중복되는 중근이 있는 경우 에는, (10)뿐 아니라 (10)을 미분한 식을 이용하면 b 개의 관계식을 구할 수 있다. 왜냐하면 (4)로부터 (9)의 분모 식 Ab(z) - zb = 0이나 분자식 (10) 모두 다항식의 형식으 로 변형될 수 있고, 다항식의 해가 중근인 경우 다항식을 미분한 식에서도 해가 되기 때문이다(자세한 내용은 Chaudhry and Templeton[5]를 참조한다.) (9)의 미지수를 구하기 위해서는 분모식의 근을 수치적으로 구해야 하는 데, 본 연구에서는 Jenkins and Traub[14] 등에서 제시한 다항식의 복소수 근을 구하는 방법을 사용하였다. 이러 한 방법 이외에도 분자와 분모에서 공통인수를 소거하는 방식으로도 미지수를 구할 수 있다. 관련된 내용은 Kim et al.[18]을 참조한다.

    취합검사 종료 직후 시점의 검체 수 분포 πi를 구하면, 취합검사의 크기(한 취합검사에 포함되는 검체 수) Y에 대한 분포와 평균을 다음처럼 구할 수 있다.

    Pr { Y = y } = { i = 0 a π i , y = a π y , a + 1 y b 1 i = b π i , y = b .
    (12)

    E [ Y ] = a i = 0 a 1 π i + y = a b 1 y π y + b i = b π i .
    (13)

    단, (13)의 무한 급수 i = b π i 는 확률의 정규 조건에서 i = b π i = 1 i = 0 b 1 π i 로 구할 수 있다.

    다음으로 임의 시점에 시스템에 있는 검체 수 L 에 대 한 분포를 구한다. 일반적으로 배치 서비스 대기행렬의 임의시점의 고객수 분포는 부가변수법(supplementary variable technique)이나 준-마아코프 과정(semi-Markov process)을 이용하여 구한다[3, 5-7, 11, 20, 22, 23, 25, 30]. 본 연구에서는 이러한 방법 대신 발생률 평형 방정 식(rate balance equations)을 이용하여 임의 시점의 검체 수 분포를 구하도록 한다. 발생률 평형 방정식을 이용하 면 부가변수법이나 준-마아코프 과정을 사용하는 경우보 다 좀 더 간단한 절차로 임의 시점의 고객수 분포를 구 할 수 있기 때문이다. 최근 Kim의 연구[17]에서 발생률 평형 방정식으로 배치 도착과 배치 서비스를 가지는 휴 가형 대기행렬 모형의 임의시점 고객 수에 대한 연구를 수행하였는데, 본 연구에서는 Kim[17]의 방법을 사용하 여 임의 시점의 검체 수 분포를 구한다.

    임의 시점의 검체 수의 분포와 PGF를 pi = Pr{L = i} 와 P ( z ) = i = 0 p i z i 로 정의하자. 그리고 λD를 취합검사 종 료 사건의 발생률(rate)이라고 하자. 임의 시점에 시스템 내의 검채 수 L 의 상태 변화는 검체의 도착으로 검체의 수가 증가하거나 취합검사 종료로 검체의 수가 감소하는 사건에 의해서만 발생한다. 먼저 L = 0이라는 상태로 진 입하는 사건을 고려해 보자. 검체의 도착으로는 L = 0이 라는 상태로 진입하지 못한다. L = 0인 상태가 되려면 취합검사가 종료 직후 시점에 시스템 내에 검체 수가 0 인 사건이 발생하여야 한다. L = i(i ≥ 1) 상태로 진입하 는 사건은 L = i - 1 상태일 때 검체가 도착하거나, 취합 검사 종료 직후 i 개의 검체가 시스템에 남아있는 경우 발생한다. 따라서 시스템 내의 검체 수 LL = i 상태 로 들어가는 진입률(rate in)은 다음과 같다.

    { λ D π 0 , i = 0 λ p i 1 + λ D π i , i 1.

    L = i로의 상태 진입률에 zi를 곱한 후 합하여 z-변환 형태로 상태 진입률을 표현하면 다음과 같이 표현된다.

    λ z P ( z ) + λ D Π ( z ) .

    마찬가지로 시스템 내의 검체 수 LL = i 상태에서 나가는 진출률(rate out)을 고려해 보자. 검체 도착으로 L = i 상태에서 나가는 진출률은 λpi이다. 취합검사 종료 로 L = i 상태에서 나가는 진출률에 대한 분석을 조금 더 복잡하다. 현재의 취합검사 종료가 m 번째 취합검사 종료하고 하자. 그러면 (m - 1) 번째 취합검사 종료 직 후의 검체 수 Lm - 1과 m 번째 취합검사 종료 직전 시점 (취합검사를 받은 고객이 아직 떠나기 직전 시점)의 고 객 수 L m 에 대한 다음의 관계가 성립한다.

    L m = { a + A a , 0 L m 1 a 1 L m 1 + A L m 1 , a L m 1 b 1 L m 1 + A b , L m 1 b .

    그러므로 안정 상태에서 L = i 상태에서의 진출률은 다음처럼 표현할 수 있다.

    λ p i + λ D Pr { lim m L m = i } .

    그리고 Lm - 1 L m 의 관계식을 이용하여 상태 진입 률을 z-변환 형태로 표현하면 다음과 같이 표현할 수 있다.

    λ P ( z ) + λ D { i = 0 a 1 π i z a A a ( z ) + i = a b 1 π i z i A i ( z ) + i = b π i z i A b ( z ) } .

    안정상태에서 모든 상태는 진출률과 진입률이 같아야 한다. 따라서 진입률과 진출률의 z-변환을 등치시키면 다 음 관계식을 얻는다.

    λ ( 1 z ) P ( z ) = λ D { n = 0 a 1 π n ( z n z a A a ( z ) ) + n = a b 1 π n z n ( 1 A n ( z ) ) + n = b π n z n ( 1 A b ( z ) ) } .

    그런데 취합검사를 받고 시스템에서 나가는 검체 이 탈률과 시스템에 들어오는 검체 도착률은 같아야 하므로 다음이 성립한다.

    λ D E [ Y ] = λ .

    따라서 앞의 두 식을 이용하면 임의 시점의 검체 수의 PGF P (z)는 취합검사 종료 직후 시점의 검체 수 분포로 다음처럼 표현할 수 있다.

    P ( z ) = 1 E [ Y ] { i = 0 a 1 π i z i z a A a ( z ) 1 z + i = a b 1 π i z i 1 A i ( z ) 1 z + i = b π i z i 1 A b ( z ) 1 z } .
    (14)

    (14)을 미분한 후 z = 1을 대입하면, 임의 시점에 시스 템에 머무르고 있는 평균 검체 수 L 을 다음과 같이 구 할 수 있다.

    L ¯ = 1 E [ Y ] [ i = 0 a 1 π i { ( a i ) ( a + i 1 ) 2 + E [ A a ] ( a + E [ A a ( A a 1 ) ] 2 E [ A a ] ) } + i = a b 1 π i E [ A i ] ( i + E [ A i ( A i 1 ) ] 2 E [ A i ] ) + i = b π i E [ A b ] ( i + E [ A b ( A b 1 ) ] 2 E [ A b ] ) ] .
    (15)

    위의 식의 Ay의 모멘트는 다음처럼 구할 수 있다.

    E [ A y ] = λ E [ S y ] , E [ A y ( A y 1 ) ] = λ 2 E [ S y 2 ] , a y b .

    또한 (15)의 무한 급수 i = b π i i 항은 다음 관계식을 이 용하여 구할 수 있다.

    i = b π i i = lim z 1 Π ( z ) i = 0 b 1 π i i .

    4. 취합검사 시스템의 성능 지표

    취합검사 시스템의 운영에는 다음과 같은 세 가지 유 형의 비용이 발생할 수 있다.

    첫째, 검체의 검사 처리가 지연되면 검사 지연에 따른 비용이 발생한다. 검사 처리 기간이 길어지면 검사 결과 를 기다리는 검사 대상자의 불편 비용(격리 기간 등)이 증가할 뿐 아니라, 양성 진단의 지연은 밀접촉자에 대한 추적을 지연시켜 감염 전파의 위험 비용을 증가시킨다. 본 연구에서는 검체 당 검사 시스템에 단위 시간 머무르 는 동안 δL의 불편 및 전파 위험 비용이 발생한다고 가 정한다.

    둘째, 검사 시스템은 유전자 증폭 검사 1회 당 검사 비용 이 발생한다. 진단 검사를 한 번 시행할 때마다 진단 키트와 시료가 사용된다. 코로나19 팬데믹의 초창기에 진단 키트와 시료를 구하지 못해 많은 국가에서 진단 검사에 어려움이 있었고, 코로나19 바이러스에 대한 취합검사도 희소한 진단 키트와 시료를 효율적으로 사용하고자 하는 목적에서 수행 되었다[21]. 아울러 진단 검사를 한 번 수행할 때마다 유전 자 증폭 검사를 하는 장비의 감가상각과 전문 인력의 투입 에 따른 인건비도 발생하게 된다. 본 연구에서는 진단 검사 1회 당 발생하는 키트와 시약비, 검사 장비의 감가 상각비, 직접 인건비 등을 단위 진단 검사 비용이라고 하고 δM이라 고 표현하기로 한다.

    셋째, 취합검사로 인해 검사의 오진단이 늘어나 사회 적 비용이 발생할 수 있다. 양성자를 음성자로 잘못 판단 하면 위음성으로 진단된 검사 대상자를 통해 감염 전파 가 이루어져 사회적 위험 비용이 발생한다. 본 연구에서 는 양섬 검체 1개가 그룹 검사에서 위음성으로 판단되어 개별 검사를 하지 않음으로써 발생하는 사회적 위험 비 용을 δF라고 가정한다.

    임의 시점에 검사 시스템에는 평균적으로 L 개의 검 체가 머무른다. 따라서 검체가 검사 시스템에 체류하는 동안 발생하는 불편 및 위험 비용은 단위 시간 당 δLL 만 큼 발생한다.

    최소 취합검사의 크기를 a, 최대 취합검사 크기를 b라 고 하면, 단위 시간 당 취합검사는 λD = λ/E [Y ]의 발생 률로 수행되고, (7), (12), (13)에서 단위 시간당 발생하는 평균 검사 횟수 M는 다음과 같이 계산될 수 있다.

    M ¯ = λ D y = a b Pr { Y = y } E [ M y ] = λ { 1 + a γ a i = 0 a 1 π i + y = a b 1 y γ y π y + b γ b i = b π i } E [ Y ] .

    단,

    γ y = k = 1 y ( y k ) p k ( 1 p ) y k γ y / k , a y b .

    따라서 단위 시간 당 발생하는 진단 검사 비용은 δMM 이 된다.

    마찬가지로 단위 시간 당 λD로 그룹 검사가 이루어지 므로, (8)과 (12)에서 단위 시간 동안 그룹 검사로 위음성 으로 판정되는 양성 검체 수의 평균 F 는 다음과 같이 계산될 수 있다.

    F ¯ = λ D y = a b Pr { Y = y } E [ F y ] = λ { ( a p γ a ) i = 0 a 1 π i + y = a b 1 ( y p γ y ) π y + ( b p γ b ) i = b π i } E [ Y ] .

    단,

    γ y = k = 1 y k ( y k ) p k ( 1 p ) y k γ y / k , a y b .

    따라서 단위 시간 당 발생하는 취합검사로 인한 위음 성 위험 비용은 δFF 가 된다.

    그러므로 취합검사 시스템의 단위 시간 당 발생하는 총 비용 Δ는 다음과 같다.

    Δ = δ L L ¯ + δ M M ¯ + δ F F ¯ .

    따라서 최적 취합검사 크기를 결정하는 문제는 다음 과 같은 최적화 문제가 된다.

    M a x i m i z e ( a , b ) Δ s . t . 1 a b b max .

    단, bmax는 기술적 한계로 취합검사를 수행할 수 있는 최대 검체 수이다.그러므로 다음 절에서 최소 및 최대 그 룹 검사 크기 ab에 따라 총 비용 Δ와 총 비용에 영 향을 주는 세 개의 성능 지표 L , M, F 가 어떻게 변화하 는지를 수치 예제를 통해 살펴보도록 한다.

    5. 수치 예제

    <Figure 1>은 의 데이터를 사용하여 2020년 3월부터 2021년 10월까지의 대한민국과 주요 선진국의 코로나19 확진율의 변화 추이를 그린 것이다. 주요국의 코로나19 확진율은 코로나19 유행이 전면화 된 후 0.2%에서 20% 사이에서 변화하였고, 대한민국의 코로나19 확진율은 1~2%의 안정된 추이를 보이다가 최근 델타 변이의 유행 으로 5% 이상으로 치솟기도 하였다. 본 연구에서는 대한 민국의 2021년의 확진율의 주요 범위를 참조하여 코로나 19 검사의 확진율이 낮은 경우(p = 0.5%), 중간인 경우 (p = 2%), 높은 경우(p = 10%)에 대하여 취합검사 수에 따른 성능 지표의 변화를 살펴볼 것이다.

    Mutesa et al.[21]의 <Table 3>에는 33명의 양성 검체 에 대한 N-유전자와 orfab를 검출에 걸린 유전자 증폭 사이클 수 Ct에 대한 데이터가 있다. 이 데이터를 분석 하면 N-유전자의 Ct는 평균 27.0, 표준편차 5.08을, orfab 의 Ct는 평균 29.02, 표준편차 5.52이고, 두 유전자의 Ct 가 정규성을 따르는지에 대한 Shapiro-Wilk 검정을 해보 면 p-값이 각각 0.6768과 0.3075로 유의수준 5%에서 정 규성에 대한 가정을 채택한다. 그러므로 본 수치 예제에 서는 양성인 개별 검체를 검출하는데 필요한 사이클의 수 Ct가 평균 μ = 29, 표준편차 σ = 5의 정규분포를 따른 다고 가정하도록 한다. 아울러 Mutesa et al.[21]과 마찬 가지로 개별 검사에서 최대 유전자 증폭 사이클 횟수는 c 1 max = 40으로 설정하도록 한다.

    마찬가지로 Mutesa et al.[21]의 연구 에서 취합검사의 검체 수가 2배 늘어날 때마다 목표 유전자 검출에 필요 한 유전자 증폭 사이클의 증가량 ΔCt가 N-유전자의 경 우 평균 1.0, 표준편차 0.15의 정규분포를, orfab의 경우 평균 1.1, 표준편차 0.14의 정규분포를 따른다고 하였으 므로, 본 수치 예제에서는 ΔCt가 평균 Δμ = 1, 표준편 차 Δσ = 0.15의 정규분포를 따른다고 가정하도록 한다. 아울러 유전자 증폭 검사의 준비 시간은 μr = 5의 지수 분포를, 유전자 증폭 한 사이클 시간은 μc = 60의 지수분 포를 따른다고 가정한다.

    본 수치 예제에서는 그룹 검사의 위음성률이 αG = 1% 이하가 되도록 c y max 를 설정한다. <Figure 2>는 확진율 p 가 0.5, 2, 10%일 때 그룹 검사의 크기 y가 커짐에 따라 c y max 의 값이 어떻게 변하는지를 보여준다(선 그래프가 겹치는 것을 방지하기 위해서 p의 값에 따라 c y max 값에 약간의 변동을 주었다.) 그룹 검사의 크기가 커지면, 그 룹 검사의 위음성률을 줄이기 위해서 유전자 증폭 사이 클의 최대 횟수 c y max 가 증가하는 것을 볼 수 있다. 특히 확진율 p가 낮을수록 한 그룹 검사에 포함되는 양성 검 체의 수가 줄어들므로 위음성률을 줄이기 위해서 c y max 가 더 빠르게 증가하는 것을 볼 수 있다.

    <Figure 3>은 그룹 검사 크기 y에 따라 검체 당 평균 취합검사 시간 E [Sy ]/y, 평균 유전자 증폭 검사 횟수 E [My ]/y, 그룹 검사에서 위음성으로 오진단 되는 평균 검 체 수 E [Fy ]/y의 변화를 보여준다. 평균 취합검사 시간 E [Sy ]는 취합검사 당 수행하는 유전자 증폭 검사 횟수 E [My ]와 밀접한 관계를 가지고 있기 때문에 검체 당 평균 취합검사 시간 E [Sy ]/y와 평균 유전자 증폭 검사 횟수 E [My ]/y 그래프는 비슷한 양상을 보인다(<Figure 3>의 왼 쪽과 가운데 그래프 참조.) 그룹 검사 크기 y가 증가하면 처음에는 검체 당 평균 취합검사 시간과 검사 횟수 모두 감소하지만, 어느 크기 이상이 되면 그룹 검사에 양성 검 체가 포함되는 확률이 급격히 증가하여 그룹 검사뿐 아니 라 개별 검사를 수행해야 하는 경우가 많아져서 검체 당 평균 검사 시간과 검사 횟수 모두 증가하게 된다. 아울러 확진율 p가 클수록 검체 당 평균 검사 시간과 횟수가 최저 가 되는 최적 그룹 검사의 크기가 작아지고, 확진율이 낮 을수록 최적의 그룹 검사 크기가 커짐을 확인할 수 있다.

    검체 당 그룹 검사에서 위음성으로 오진단되는 평균 검체 수 E [Fy ]/y는 그룹 검사 수 y가 변함에 따라 톱니바퀴 모양 의 경향을 보인다(<Figure 3>의 오른쪽 그래프 참조.). 이러 한 경향을 보이는 이유는, <Figure 2>에서 보듯이 유전자 증폭 횟수 c y max 의 값이 y가 변함에 따라 계단식으로 변화하 기 때문이다. 그룹 검사 크기 y의 값이 어느 정도 커지면, 동일한 유전자 증폭 횟수 c y max 를 가지는 y에 대하여 E [Fy ]/y 는 증가하는 경향을 보이다가, c y max 의 값이 바뀌는 y에서 크게 하락하는 패턴을 보인다. 동일한 c y max 값을 가지는 y에 대하여 E [Fy ]/y가 증가하는 경향을 보이는 이유는 다음 과 같다. (1)에서 보듯이 코로나19 유전체가 검출되는 유전 자 증폭 사이클 수 CTY/K의 분포는 log(Y/K ) 에 비례하여 평균과 표준편차가 증가한다. 따라서 동일한 양성 검체 수 K를 포함하는 그룹 검사는 Y가 커짐에 따라 그룹 검체의 농도가 묽어져서 검출 확률이 줄어든다. 그러므로 y가 커짐 에 따라 그룹 위음성률이 커지게 되고, 이 위음성률이 어느 정도 이상이 되면 c y max 의 값이 조정되어 다시 감소하는 지그 재그 패턴을 보이게 된다.

    <Figure 4>는 개별 검체의 도착률이 λ = 0.4이고 최대 그룹 검사 크기가 b = 10, 그리고 검체 하나가 시스템에 단위 시간 동안 대기하는 동안 발생하는 불편 비용이 δL = 1, 검사 당 발생하는 진단 비용이 δM = 50, 위음성으 로 오진하여 발생하는 위험 비용이 δF = 1000일 때, 최소 그룹 검사 크기 a에 따라 취합검사 시스템의 성능 지표 가 어떻게 변화하는지를 보여준다. <Figure 4>의 좌상단 의 그래프는 시스템에 머무는 평균 검체 수 L 의 변화를 보여준다. 최소 그룹 검사 크기 a가 증가하면 확진율 p 에 무관하게 증가하는 경향을 보여준다. 이는 최소 그룹 검사 크기가 커지면 최소 그룹 검사 크기가 될 때까지 대기해야 하는 시간이 커지기 때문에 발생하는 현상이 다. <Figure 4>의 우상단의 그래프는 단위 시간 당 시스 템에서 수행해야 하는 평균 유전자 증폭 검사 수 M의 변화를 보여준다. 최소 그룹 검사 크기 a가 증가하면 처 음에는 그룹 검사의 효과로 평균 검사 수가 감소한다. 그 러나 확진율이 높을 때에는 최소 검사 수가 어느 값 이 상 커지면 그룹 검사 후에 개별 검사를 수행하는 경우가 많아져 오히려 평균 검사 수가 증가하는 것을 보여준다 (<Figure 3>의 가운데 그래프 참조.) <Figure 4>의 좌하 단의 그래프는 단위 시간 당 그룹 검사에 의해 위음성으 로 판정되어 검역 체계를 벗어나는 평균 양성 검체 수 F 의 변화를 보여준다. 그룹 검사에서 위음성으로 오진단 되는 평균 양성 검체 수는 a의 변화에 따라 지그재그 형 식으로 변화하는 것을 볼 수 있다. 이러한 현상이 나타나 는 이유는, <Figure 3>의 오른쪽 그래프에 대한 설명에서 언급한 것처럼, 그룹 검사는 그룹 검사 크기가 커짐에 따 라 그룹 검사의 위음성률이 커지게 되고, 이 위음성률이 어느 정도 이상이 되면 c y max 의 값이 조정되어 다시 감소 하는 지그재그 패턴을 보이기 때문이다. 마지막으로 <Figure 4>의 우하단의 그래프는 취합검사 시스템의 총 비용의 변화를 보여준다. 최소 그룹 검사 크기가 커지면 대기 시간이 길어져 비용이 증가하지만, 한 번의 그룹 검 사로 여러 번의 개별 검사를 대체할 수 있으므로 검사 비용은 감소하게 된다. 총비용 함수는 이러한 두 요인의 상쇄 효과에 의해 아래로 오목한 곡선의 형태를 보인다. 그렇기 때문에 총비용이 최소가 되는 지점에서 최적인 최소 그룹 검사 크기 a*를 설정할 수 있다. 확진율 p가 커지면 그룹 검사의 효율이 나빠지기 때문에 최적인 최 소 그룹 검사 크기가 작아지는 것을 볼 수 있다.

    <Figure 5>는 개별 검체의 도착률과 비용 요인이 동일 하고 최소 그룹 검사 크기가 a = 2일 때, 최대 그룹 검사 크기 b가 변하면 취합검사 시스템의 성능 지표가 어떻게 변화하는지를 보여준다. 모든 성능 지표에서 최소 그룹 검사 크기 a의 그래프인 <Figure 4>에 비해 변화가 크지 않음을 볼 수 있다. 이는 시스템의 처리 용량보다 검체 도착률이 낮으면 최대 그룹 검사 크기에 이르기 전에 서 비스가 이루어져서 최대 그룹 검사 크기는 성능에 크게 영향을 미치지 않기 때문이다.

    검사 요청의 도착률 λ가 커지면 최대 그룹 검사 크기 b가 시스템 성능 지표에 영향을 크게 주는 지를 살펴보기 위하여 λ의 값을 0.4, 0.8, 1.2로 바꾸어 가며 b에 따른 시 스템의 성능 지표의 변화를 살펴보자. <Figure 6>은 중간 정도의 확진율 p = 2%일 때 검체 도착률 λ의 세 가지 값 에 대해 최대 그룹 검사 크기 b에 따른 성능 지표의 변화 를 시각화한 것이다. <Figure 6>의 좌상단의 그래프에서 보듯이 검체 도착률이 큰 경우 최대 그룹 검사 크기가 커 지면 검체의 평균 대기 시간을 감소시기는 효과가 나타나 지만, 최대 그룹 검사 크기가 어느 정도 늘어나면 이러한 효과는 정체되는 것을 볼 수 있다. 그 이유는 앞서 <Figure 5>에 대한 설명에서도 언급하였듯이, 최대 그룹 검사 크기 가 어느 정도 커져 시스템의 처리 용량이 도착하는 검체의 양보다 커지면 최대 그룹 검사에 도달하기 전에 처리되는 경우가 많아져서 최대 그룹 검사 크기는 시스템 성능에 큰 영향을 미치지 못하기 때문이다.

    지금까지의 수치 예제의 분석을 통해 다음과 같은 결 론을 얻을 수 있다. 첫째, 취합검사 시스템의 용량이 도 착하는 검사 요구량을 충분히 감당할 수 있는 경우, 취합 검사 시스템의 성능은 최대 그룹 검사 크기보다는 최소 그룹 검사 크기에 영향을 더 크게 받는 경향을 보인다. 둘째, 검사의 확진율이 높으면 최소 검사 크기를 작게 조 정하고, 확진율이 낮으면 최소 검사 크기를 크게 조정하 는 것이 검체 대기에 따른 불편 비용, 검사 비용, 위음성 진단에 따른 위험 비용을 모두 고려한 총 비용을 낮출 가능성이 커진다. 셋째, 검사 비용과 위음성률만 고려할 때보다 검체의 시스템 대기 시간을 고려하면 최적인 최 소 취합검사 크기가 더 작아지는 경향을 보인다. 왜냐하 면 검체의 시스템 대기 시간은 최소 취합검사 크기가 커 지면 지속적으로 커지기 때문에(<Figure 4> 좌상단 그래 프 참조), 총비용이 최소가 되는 최소 취합검사 크기가 더 낮은 영역에서 형성되기 때문이다.

    6. 결 론

    도프만이 취합검사법을 제안한 후, 취합검사의 성능을 분석하는 많은 연구들이 있었다. 그러나 대부분의 연구 들은 검사 요구의 도착 과정의 확률적 속성과 검사 시간 의 확률적 속성을 고려하지 않고, 확정적인 상황을 가정 하여 취합검사 시스템을 분석하였다. 검사 요구 도착 과 정과 검사 시간의 확률적 속성을 고려한 연구로는 시뮬 레이션 기법을 사용한 연구와 대기행렬을 이용한 연구가 있다. 시뮬레이션을 이용한 연구는 현실의 취합검사 시 스템을 유사하게 모사하여 연구할 수 있는 장점이 있지 만, 다양한 시나리오를 고려해야 할 때 시뮬레이션 시간 이 길어지는 단점이 있다. 반면 대기행렬을 이용한 연구 는 분석적 방법으로 시스템을 탐색할 수 있는 장점이 있 지만, 지금까지 대기행렬을 이용한 취합검사 시스템에 대한 연구는 검사가 항상 정확하여 위음성 진단이 없다 고 가정하였다. 그러나 코로나19 검사처럼 위음성 가능 성이 있는 경우에는, 취합검사 시스템의 성능 분석에 위 음성 가능성을 고려할 필요가 있다. 그러므로 본 연구에 서는 배치 서비스 대기행렬을 이용하여 위음성률이 있는 코로나19 취합검사 시스템의 성능을 분석하는 방법을 제 시하였다. 아울러 코로나19 취합검사에 대한 기존 연구 를 반영하여, 취합검사 시간, 검사 횟수, 위음성으로 진 단되는 양성 검체 수의 분포를 모형화였고, 수치 예제를 통해 확진율과 취합검사의 크기에 따라 취합검사 시스템 의 성능이 어떻게 변하는지 살펴보았다. 이를 통해 시스 템의 확률적 속성 때문에 발생하는 검체의 처리 시간을 고려할 때, 최대 취합검사 크기보다는 최소 취합검사 크 기가 시스템 성능에 더 큰 영향을 미치는 것을 살펴보았 고, 최적 최소 취합검사 크기는 검체의 처리 시간을 고려 하는 경우가 이를 고려하지 않을 때보다 작아지는 경향 이 있다는 것을 예시하였다. 이는 코로나19처럼 신속한 검사 결과가 필요한 경우 취합 검사의 개시 크기를 더 작게 설정하는 것이 필요하다는 것을 시사한다.

    마지막으로 본 연구의 향후 발전방향을 제시하며 논문을 마무리하고자 한다. 첫째, 코로나19 유전자 증폭 검사의 위양성률은 0%이기 때문에, 본 연구에서는 위음성률만 있 는 취합검사 시스템에 대해서 분석하였다. 그러나 다른 검 사 시스템은 위음성뿐 아니라 위양성 판단이 있을 수 있기 때문에 위양성과 위음성이 모두 있는 취합검사 시스템을 대기행렬로 분석하는 연구로 확장될 수 있을 것이다. 둘째, 본 연구에서는 코로나19 취합 검사 시스템을 예로 분석하였 지만, 검사 요구의 도착과 검사 시간이 확률적 속성이 있고 위음성이나 위양성 진단이 있는 다른 취합 진단 시스템에도 연구 방법이 확대 적용 가능할 것으로 기대된다. 셋째, 본 연구의 수치 예제에서는 유전자 증폭 검사의 사이클 수의 분포 등은 기존 연구를 참조하여 현실적인 수치로 분석을 하였지만, 검사 요구의 도착률, 비용 요소 등은 임의의 수치 가 사용되었다. 향후 실제 검사 시스템의 데이터를 측정하 여 분석한다면 좀 더 현실적인 연구가 될 수 있을 것이라 기대된다.

    Figure

    JKISE-44-4-154_F1.gif

    The Positive Rates of Daily COVID-19 Tests for South Korea and Other Countries

    JKISE-44-4-154_F2.gif

    Cymax for Various Values of p

    JKISE-44-4-154_F3.gif

    E [Sy ]/y, E [ My ]/y, and E [Fy]/y over the group test size y for various values of p

    JKISE-44-4-154_F4.gif

    Performance Measures over the Minimum Group Test Size a for Various Values of p

    JKISE-44-4-154_F5.gif

    Performance Measures over the Maximum Group Test Size b for Various Values of p

    JKISE-44-4-154_F6.gif

    Performance Measures over the Maximum Group Test Size b for Various Values of λ

    Table

    Reference

    1. Abdalhamid, B., Bilder, C.R., McCutchen, E.L., Hinrichs, S.H., Koepsell, S.A., and Iwen, P.C., Assessment of specimen pooling to conserve SARS CoV-2 testing resources, American Journal of Clinical Pathology, 2020, Vol. 153, No. 6, pp. 715-718.
    2. Abolnikov, L. and Dukhovny, A., Optimization in HIV screening problems, Journal of Applied Mathematics and Stochastic Analysis, 2003, Vol. 16, No. 4, pp. 361-374.
    3. Bailey, N.T., On queueing processes with bulk service, Journal of the Royal Statistical Society: Series B (Methodological), 1954, Vol. 16, No. 1, pp. 80-87.
    4. Bar-Lev, S.K., Parlar, M., Perry, D., Stadje, W., and Van der Duyn Schouten, F.A., Applications of bulk queues to group testing models with incomplete identification, European Journal of Operational Research, 2007, Vol. 183, No. 1, pp. 226-237.
    5. Chaudhry, M. and Templeton, J.G., First Course in Bulk Queues, John Wiley & Sons, 1983.
    6. Claeys, D., Steyaert, B., Walraevens, J., Laevens, K., and Bruneel, H., Tail probabilities of the delay in a batch-service queueing model with batch-size dependent service times and a timer mechanism, Computers & Operations Research, 2013, Vol. 40, No. 5, pp. 1497-1505.
    7. Claeys, D., Walraevens, J., Laevens, K., and Bruneel, H., A queueing model for general group screening policies and dynamic item arrivals, European Journal of Operational Research, 2010, Vol. 207, No. 2, pp. 827-835.
    8. Claeys, D., Walraevens, J., Laevens, K., and Bruneel, H., Analysis of threshold-based batch-service queueing systems with batch arrivals and general service times, Performance Evaluation, 2011, Vol. 68, No. 6, pp. 528-549.
    9. Deckert, A., Bärnighausen, T., and Kyei, N.N., Simulation of pooled-sample analysis strategies for COVID-19 mass testing, Bulletin of the World Health Organization, 2020, Vol. 98, No. 9, pp. 590.
    10. Dorfman, R., The detection of defective members of large populations, The Annals of Mathematical Statistics, 1943, Vol. 14, No. 4, pp. 436-440.
    11. Downton, F., Waiting time in bulk service queues, Journal of the Royal Statistical Society: Series B (Methodological), 1955, Vol. 17, No. 2, pp. 256-261.
    12. Du, D., Hwang, F.K. and Hwang, F., Combinatorial group testing and its applications, World Scientific, 2000.
    13. Garg, J., Singh, V., Pandey, P., Verma, A., Sen, M., Das, A., and Agarwal, J., Evaluation of sample pooling for diagnosis of COVID-19 by real time-PCR: A resource-saving combat strategy, Journal of Medical Virology, 2021, Vol. 93, No. 3, pp. 1526-1531.
    14. Jenkins, M.A. and Traub, J.F., Algorithm 419: zeros of a complex polynomial, Communications of the ACM, 1972, Vol. 15, No. 2, pp. 97-99
    15. Kim, J.-H., Comparative study of target genes and protocols by country for detection of SARS-CoV-2 based on polymerase chain reaction (PCR), Journal of the Korea Contents Association, 2021, Vol. 21, No. 1, pp. 465-474.
    16. Kim, J. and Lee, T., A study on the application and optimization of pooled tests for the diagnosis of COVID-19, Proceedings of the 2021 KIIE (Korean Institute of Industrial Engineers) Spring Conference, 2021, pp. 4516-4554.
    17. Kim, K., A relation between queue-length distributions during server vacations in queues with batch arrivals, batch services, or multiclass arrivals: An extension of burke’s theorem, Indian Journal of Science and Technology, 2015, Vol. 8, No. 25.
    18. Kim, N.K., Chaudhry, M.L., Yoon, B.K. and Kim, K., Inverting generating functions with increased numerical precision: A computational experience, Journal of Systems Science and Systems Engineering, 2011, Vol. 20, pp. 475-494.
    19. Macula, A.J., Probabilistic nonadaptive group testing in the presence of errors and DNA library screening, Annals of Combinatorics, 1999, Vol. 3, No. 1, pp. 61-69.
    20. Medhi, J., Waiting time distribution in a poisson queue with a general bulk service rule, Management Science, 1975, Vol. 21, No. 7, pp. 777-782.
    21. Mutesa, L., Ndishimye, P., Butera, Y., Souopgui, J., Uwineza, A., Rutayisire, R., Ndoricimpaye, E.L., Musoni, E., Rujeni, N., Nyatanyi, T. and others, A pooled testing strategy for identifying SARS-CoV-2 at low prevalence, Nature, 2021, Vol. 589, No. 7841, pp. 276-280.
    22. Nandy, N. and Pradhan, S., On the joint distribution of an infinite-buffer discrete-time batch-size-dependent service queue with single and multiple vacations, Quality Technology & Quantitative Management, 2021, pp. 1-36.
    23. Neuts, M.F., A general class of bulk queues with poisson input, The Annals of Mathematical Statistics, 1967, Vol. 38, No. 3, pp. 759-770.
    24. Ngo, H.Q. and Du, D.-Z., A survey on combinatorial group testing algorithms with applications to DNA library screening, Discrete Mathematical Problems with Medical Applications, 2000, Vol. 55, pp. 171-182.
    25. Pradhan, S. and Gupta, U., Analysis of an infinite-buffer batch-size-dependent service queue with markovian arrival process, Annals of Operations Research, 2019, Vol. 277, No. 2, pp. 161-196.
    26. Song, G.S., Lee, Y.-R., Kim, S., Kim, W., Choi, J., Yoo, D., Yoo, J., Jang, K.-T., Lee, J. and Jun, J.H., Laboratory diagnosis of coronavirus disease 19 (COVID-19) in korea: Current status, limitation, and challenges, Korean J Clin Lab Sci, 2020, Vol. 52, No. 3, pp. 284-295.
    27. Wein, L.M. and Zenios, S.A., Pooled testing for HIV screening: Capturing the dilution effect, Operations Research, 1996, Vol. 44, No. 4, pp. 543-569.
    28. Yelin, I., Aharony, N., Tamar, E.S., Argoetti, A., Messer, E., Berenbaum, D., Shafran, E., Kuzli, A., Gandali, N., Shkedi, O. and others, Evaluation of COVID-19 RT-qPCR test in multi sample pools, Clinical Infectious Diseases, 2020, Vol. 71, No. 16, pp. 2073-2078.
    29. Yong, D., Laboratory diagnosis of 2019 novel coronavirus, Korean J healthc assoc Infect Control Prev, 2020, Vol. 25, No. 1, pp. 63-65.
    30. Zhao, Y.Q. and Campbell, L.L., Equilibrium probability calculations for a discrete-time bulk queue model, Queueing Systems, 1996, Vol. 22, No. 1, pp. 189-198.